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14.  ABSTRACT 

This  report  describes  the  development  and  evaluation  of  a  cepstral  analysis,  modeling  and  inversion  program  for  the  characterization  of 
underwater  explosions  recorded  by  seismic  stations.  The  algorithm  matches  synthetic  against  observed  cepstrums  for  suspected 
underwater  blasts.  The  observed  cepstrums  are  computed  for  regional  phases  (e.g.,  Pn,  Pg,  Sn,  and  Lg)  by  taking  the  logarithm  of  the 
trend-corrected  spectrums  for  each  phase,  stacked  across  an  regional  array  if  array  data  is  available,  and  then  taking  the  inverse  Fourier 
transform  of  log  amplitude  spectrum.  The  signed  cepstrums  for  regional  phases  from  underwater  explosions  have  negative  peaks  caused 
by  the  reflections  of  the  acoustic  wave  from  the  surface  and  positive  peaks  from  the  bubble  pulse.  The  depth  and  yield  of  the  underwater 
blast  are  determined  by  finding  a  synthetic  cepstrum  that  most  closely  matches  the  observed  cepstrum.  The  best  matching  synthetic 
cepstrum  is  the  one  with  the  highest  match,  either  by  cross  correlation  coefficient,  LI  Norm,  or  L2  Norm,  with  the  observed  cepstrum.. 
This  algorithm  has  been  tested  on  the  November  1999  DTRA  sponsored  calibration  blasts  in  the  Dead  Sea,  and  the  method  gives 
explosion  yields  consistent  with  the  known  yields  and  depths  of  the  events.  Results  are  also  given  for  analysis  of  two  events  in  the 
Barents  Sea,  an  explosion  near  Murmansk  and  the  Kursk  submarine  event.  We  conclude  that  the  Kursk  event  was  an  underwater 
explosion  with  a  yield  of  about  4300  kg,  or  4.73  tons,  at  a  depth  of  90  m. 
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SUMMARY 


This  study  has  investigated  using  seismic  data  to  characterize  underwater  explosions  and 
discriminate  underwater  blasts  from  earthquakes.  Seismic  recordings  of  underwater  blasts 
can  be  used  to  characterize  the  source  itself,  based  on  identification  of  bubble  pulses,  and 
the  in-situ  water  depths  from  mode-converted  acoustic  signals  reflecting  in  the  vicinity  of 
the  source.  Thus,  regional  phases  from  underwater  events,  such  as  Pn,  Sn,  and  Lg,  may 
provide  more  information  about  the  source  and  in-situ  source  characteristics  than  later 
arriving  T  phases,  usually  recorded  hydroacoustically,  since  the  latter  signals  may  be 
strongly  distorted  by  heterogeneities  in  the  water-column  propagation  path.  Bubble  pulses 
and  near-source  water-column  reflections  cause  time-independent  scalloping  of  the  spec¬ 
tra  of  regional  phases  that  can  be  analyzed  to  infer  depth  and  yield  of  the  explosion. 

An  inversion  algorithm  has  been  prototyped  in  Matlab©  to  characterize  the  regional- 
phase  cepstra  from  underwater  explosions  recorded  at  seismic  stations.  The  algorithm 
matches  synthetic  against  observed  cepstrums  for  suspected  underwater  blasts.  The  ob¬ 
served  cepstrums  are  computed  for  regional  phases  (e.g.,  Pn,  Pg,  Sn,  and  Lg)  by  taking 
the  logarithm  of  the  trend-corrected  spectrums  for  each  phase,  stacked  across  an  regional 
array  if  array  data  is  available,  and  then  taking  the  inverse  Fourier  transform  of  log  ampli¬ 
tude  spectrum.  The  result  is  the  so-called  signed  cepstrum.  Each  of  the  individual  phase 
cepstrums  is  then  stacked  to  produce  a  composite  cepstrum  for  the  event.  The  signed  cep¬ 
strums  for  regional  phases  from  underwater  explosions  have  negative  peaks  caused  by  the 
reflections  of  the  acoustic  wave  from  the  surface  and  positive  peaks  from  the  bubble 
pulse.  Undersea  earthquakes  may  only  have  the  negative  water-column  reflection  peaks  or 
no  peaks  at  all.  Cepstrums  of  the  background  noise  before  Pn  are  also  computed  in  order 
to  check  if  the  peaks  are  caused  by  noise  or  processing  artifact.  Synthetic  cepstrums  are 
computed  for  assumed  minimum-phase  wavelets  for  an  explosion  of  specified  yield  and 
depth  in  the  water,  which  can  be  constrained  by  known  bathymetry.  The  depth  and  yield 
of  the  underwater  blast  are  determined  by  finding  a  synthetic  cepstrum  that  most  closely 
matches  the  observed  cepstrum.  The  best  matching  synthetic  cepstrum  is  the  one  with  the 


highest  match,  either  by  cross  correlation  coefficient,  LI  Norm,  or  L2  Norm,  with  the  ob¬ 
served  cepstrum.  The  best  matching  cepstrum  can  be  found  by  either  exhaustive  and  op¬ 
timal  search  algorithms,  including  the  downhill  simplex  method.  This  algorithm  can  iden¬ 
tify  undersea  events  as  either  blast  or  earthquake  based  on  the  degree  of  correlation  be¬ 
tween  a  best  fitting  underwater  explosion  model  cepstrum  and  the  observed  cepstrum. 

The  algorithm  has  been  tested  on  offshore  events  that  are  suspected  to  be  underwater 
chemical  blasts.  Offshore  events  have  been  collected  by  a  survey  of  the  Reviewed  Event 
Bulletins  (REB)  of  the  Prototype  National  Data  Center  (PIDC).  The  main  focus  initially 
was  on  events  around  Scandinavia,  in  the  Norwegian,  North,  and  Baltic  Seas,  and  the 
Gulf  of  Bothnia,  recorded  at  one  or  more  of  the  Scandinavian  arrays.  In  this  report,  we 
test  the  cepstral  inversion  algorithm  on  the  calibrated  Dead  Sea  explosion  data,  and  apply 
the  algorithm  to  events  in  the  Barents  Sea,  a  Russian  underwater  blast  near  Murmansk 
and  an  event  that  coincided  with  the  Kursk  submarine  disaster.  The  three  calibrated  Dead 
Sea  explosions  had  actual  yields  of  500,  2000  and  5000  kg.  Our  inversion  program  picked 
estimated  yields  of  600,  1500,  and  3100  kg,  respectively.  For  the  Kursk  event,  we  deter¬ 
mined  that  the  second  large  event  was  an  underwater  explosion.  Our  inversion  program 
picked  a  model  that  has  a  yield  of  4300  kg,  or  about  4.73  tons,  and  a  depth  of  90  m.  This 
supports  news  reports  of  a  very  large  explosions  being  related  to  the  Kursk  submarine 
tragedy. 

Although  the  cepstral-matching  algorithm  seems  to  work  well  in  characterizing  underwa¬ 
ter  blasts,  there  have  been  a  couple  of  observations  that  have  been  difficult  to  explain. 
First,  our  estimated  yields  from  bubble  pulses  to  not  correlate  strongly  with  the  local 
magnitude  estimates  for  blasts  in  different  regions.  We  noted  some  correlation  for  several 
events  in  about  the  same  location.  The  lack  of  correlation  with  extreme  values  of  yield 
might  be  due  to  confusion  in  identifying  the  true  first  bubble  pulse.  Second,  the  models 
that  fit  the  observed  cepstra  the  best  must  have  very  small  surface  reflection  coefficients, 
which  are  on  the  order  of -0.2  to  -0.4,  much  less  than  the  expected  value  of— 1.  This  re¬ 
sult  may  be  due  to  poor  resolution  of  the  peaks  at  low  queffency. 

However,  we  have  found  that  the  method  overall  is  robust  and  provides  a  means  for  char- 


iii 


acterizing  and  screening  undemater  explosions.  A  uses  manual  for  the  algorithm  and  user 
interfaces,  prototyped  in  Matlab©,  is  included  in  the  Appendix. 


CONVERSION  TABLE 

Conversion  Factors  for  U.S.  Customary  to  metric  (SI)  units  of  measurement. 


MULTIPLY - ►BY  - - - ►  TO  GET 

TO  GET  ◄ - BY  ◄ -  DIVIDE 


angstrom 

1.000  000  x  E  -10 

metexs  (m) 

atmosphere  (normal) 

1.013  25  x  E  +2 

kilo  pascal  (kPa) 

bar 

1.000  000  x  E  +2 

kilo  pascal  (kPa) 

bam 

1.000  000  X  E  -28 

meter2  (m2) 

British  thermal  unit  (thermochemical) 

1.054  350  x  E  +3 

joule  (J) 

calorie  ( thermochemical ) 

4.184  000 

joule  (J) 

cal  ( thermochemical /  cm2 ) 

4.184  000  x  E  -2 

mega  joule/m2  (MJ/m2) 

curie 

3.700  000  X  E  +1 

*giga  bacquerel  (GBq) 

degree  (angle) 

1.745  329  x  E  -2 

radian  (rad) 

degree  Fahrenheit 

tjc  =  (t°f  +  459. 67)  /l.  8 

degree  kelvin  (K) 

electron  volt 

1.602  19  x  E  -19 

joule  (J) 

erg 

1.000  000  x  E  -7 

joule  (J) 

erg/ second 

1.000  000  x  E  -7 

watt  (W) 

foot 

3.048  000  x  E  -1 

meter  (m) 

foot  -pound-  force 

1.355  818 

joule  (J) 

gallon  (U.S.  liquid) 

3.785  412  X  E  -3 

meter3  (m?) 

inch 

2.540  000  x  E  -2 

meter  (m) 

jerk 

1.000  000  x  E  +9 

j  oule  ( J) 

joule/kilcgram  (J/kg)  radiation  dose 

absorbed 

1.000  000 

Gray  (Gy) 

kilotons 

4.183 

tera joules 

kip  (1000  Ibf) 

4.448  222  X  E  +3 

newton  (N) 

kip/ inch2  (ksi) 

6.894  757  x  E  +3 

kilo  pascal  (kPa) 

ktap 

1.000  000  x  E  +2 

newton  -  second/m2  (N  -  s  /m2 ) 

micron 

1.000  000  x  E  -6 

meter  (m) 

mil 

2.540  000  X  E  -5 

meter  (m) 

mile  ( intemat ional ) 

1.609  344  x  E  +3 

meter  (m) 

ounce 

2.834  952  X  E  -2 

kilogram  (kg) 

pound- force  (lbs  avoirdupois) 

4.448  222 

newton  (N) 

pound- force  inch 

1.129  848  x  E  -1 

newton -meter  (N-m) 

pound- force /  inch 

1.751  268  x  E  +2 

newton/ me  ter  (N/m) 

pound - force/ foot2 

4.788  026  x  E  -2 

kilo  pascal  (kPa) 

pound- force/ inch2  (psi) 

6.894  757  1 

kilo  pascal  (kPa) 

pound -mass  (Ibm  avoirdupois)  1 

4.535  924  x  E  -1 

kilogram  (kg) 

pound-mass -foot2  (moment  of  inertia) 

4.214  011  x  E  -2 

kilogram-meter2  (kg-m2) 

pound -mass/  foot3 

1.601  846  X  E  +1 

kilogram-meter3  (kg/m3) 

rad  (radiation  dose  absorbed) 

1.000  00 0  x  E  -2 

**Gray  (Gy) 

roentgen 

2.579  760  X  E  -4 

coulcmb/kilogram  (C/kg) 

shake 

1.000  000  x  E  -8 

second  (s) 

slug 

1.459  390  x  E  +1 

kilogram  (kg) 

torr  (mm  Hg,  0°  C) 

1.333  22  x  E  -1 

kilo  pascal  (kPa) 

*The  bacquerel  (Bq)  is  the  SI  unit  of  radioactivity;  1  Bq  =  1  event/s. 
**The  Gray  (GY)  is  the  SI  unit  of  absorbed  radiation. 
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SECTION  1 
INTRODUCTION 


The  problem  of  the  identification  of  underwater  blasts  has  gained  increased  interest  re¬ 
cently  in  the  context  of  the  monitoring  of  the  Comprehensive  Nuclear  Test-Ban  Treaty 
(CTBT),  which  was  opened  for  signature  by  the  United  Nations  on  24  September  1996. 
Annex  1  to  the  Protocol  of  the  CTBT  calls  for  the  installation  of  an  International  Moni¬ 
toring  System  (IMS)  including  six  hydroacoustic  stations  and  five  so-called  “T -phase” 
stations.  T-phase  stations  are  seismic  sensors  located  near  the  coast  that  can  detect  hy¬ 
droacoustic  phases  converted  to  seismic  phases  at  the  coast.  Thus,  only  1 1  stations  will  be 
available  specifically  for  monitoring  underwater  events.  If  an  explosion  occurs  in  the 
ocean,  but  near  the  coast  outside  of  the  SOFAR  channel,  long-range  propagation  of  hy¬ 
droacoustic  signals  may  be  blocked,  and  there  is  a  possibility  that  the  events  may  not  be 
easily  detected  by  the  IMS  hydroacoustic  and  T-phase  assets  directed  toward  the  under¬ 
water  explosions.  Because  of  the  relatively  larger  number  of  seismic  stations,  170  primary 
plus  auxiliary  stations,  called  for  in  the  CTBT  Protocol  for  the  IMS,  near-coast  seismic 
stations  may  have  a  better  chance  of  detecting  and  characterizing  underwater  events  on 
the  continental  slopes,  outside  of  the  SOFAR  channel,  or  in  confined  seas.  Moreover, 
early-arriving  seismic  signals,  such  as  Pn,  Pg,  Sn ,  and  Lg,  produced  by  mode  conversion 
of  acoustic  waves  in  the  water  in  the  vicinity  of  the  source,  may  carry  more  useful  infor¬ 
mation  about  in-situ  source  conditions  than  later  arriving  T  phases  that  may  be  affected 
by  propagation  path  effects  in  the  oceanic  water  column. 

Recent  studies  of  underwater  acoustics  data  for  CTBT  monitoring  have  focused  primarily 
on  analysis  of  T  phase  signals  from  underwater  events.  T  phases  from  undersea  earth¬ 
quakes  are  usually  the  largest  signals  recorded  on  hydrophones  (e.g.,  deGroot-Hedlin  and 
Orcutt,  1997).  Gitterman  et  al  (1994)  have  analyzed  seismic  recordings  of  explosions  in 
the  Dead  Sea,  and  describe  a  low-frequency  spectral  analysis  method  (LFSM)  for  dis¬ 
crimination  of  underwater  explosions  from  earthquakes.  Laney  et  al,  1999  have  used  cep- 
stral  and  cross  correlation  analysis  to  characterize  bubble  pulse  signals  from  calibration 
explosions  recorded  at  hydrophone  sensors.  Calculations  of  first  bubble  pulse  period  at 
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the  P1DC,  using  cross  correlation  and  cepstral  analysis,  have  agreed  with  theoretical  val¬ 
ues  from  Urick  (1983)  for  some  calibration  explosions  (Ronan  Le  Bras,  personal  commu¬ 
nication). 

Baumgardt  and  Der  (1998)  and  Baumgardt  (1999)  published  an  extensive  study  of  nu¬ 
merous  examples  of  underwater  explosions  recorded  at  seismic  stations  and  how  they  can 
be  characterized  by  spectral  and  cepstral  analysis.  A  simple  model  for  underwater  explo¬ 
sions  was  developed  and  synthetic  cepstra  were  produced  that  reproduced  most  of  the  es¬ 
sential  features  of  observed  underwater  explosion  cepstra.  The  main  features  were  bubble 
pulses,  which  produce  positive  cepstral  peaks,  and  the  surface  reflections  that  produce 
strong  negative  cepstral  peaks.  The  timing  and  relative  amplitudes  of  these  cepstral  peaks 
provide  useful  constraints  on  the  depth  and  yield  of  underwater  explosions. 

In  the  study  described  in  this  report,  we  have  extended  our  original  work  in  Baumgardt 
and  Der  (1998)  to  directly  match  model  cepstra  to  observed  cepstra  with  the  intent  of 
identifying  and  characterizing  underwater  explosions.  Our  method  of  “cepstral  inver¬ 
sion”  involves  both  exhaustive  and  optimal  search  algorithms  that  find  models  that  fit  the 
observed  cepstra  by  maximizing  a  fit  parameter  between  the  observed  and  theoretical 
cepstra.  In  essence,  modeling  cepstra  involves  modeling  the  correlation  structure  of  the 
different  pulses  produced  in  the  water  by  underwater  explosions. 

This  report  describes  the  method  of  using  cepstral  analysis  to  characterize  spectral  modu¬ 
lations  of  underwater  explosions,  and  examples  of  its  application  to  a  presumed  underwa¬ 
ter  explosion  off  the  coast  of  Norway  near  Tromso,  the  Dead  Sea  calibration  explosions,  a 
Russian  underwater  explosion  near  Murmansk,  and  the  recent  Kursk  Submarine  event  in 
the  Barents  Sea.  The  theory  for  modeling  cepstrums  for  underwater  explosions,  originally 
described  by  Baumgardt  and  Der  (1998),  and  our  method  fitting  a  best-fitting  model  to 
observed  cepstra,  first  described  by  Baumgardt  (1999),  will  be  reviewed.  The  last  section 
will  include  summary  and  conclusions. 
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SECTION  2 

UNDERWATER  EXPLOSION  MODELING  APPROACH 


The  most  striking  feature  of  the  observations  of  underwater  explosions  is  the  occurrence 
of  strong  modulations  or  scalloping  in  the  spectra  of  each  regional  phase  that  produces 
sharp  cepstral  peaks.  Similar  features  have  been  observed  in  hydroacoustic  data  from  un¬ 
derwater  explosions  caused  by  the  complex  interference  of  bubble  pulses  and  water- 
column  bounces.  In  this  section,  we  present  a  model  and  synthetic  cepstra  that  explain  the 
seismic  observations,  and  estimate  approximate  ranges  of  sizes  and  water  depths  of  the 
charges  that  could  have  produced  the  observed  cepstra. 

An  underwater  explosion  produces  the  seismic  signals  by  the  conversion  at  the  water- 
bottom  interface  of  the  acoustic  waves  in  the  water  into  each  of  the  regional  seismic 
phases  that  propagate  through  the  earth  to  the  stations.  The  strength  of  this  conversion 
depends  on  the  transfer  function  at  the  water  bottom,  and  this  function  will  affect  the  rela¬ 
tive  strengths  of  the  different  regional  phases,  Pn,  Pg,  Sn ,  and  Lg,  which  in  turn,  controls 
the  size  of  the  P/S  amplitude  ratios.  However,  the  spectral  characteristics  that  are  com¬ 
mon  to  each  of  the  phases  are  controlled  strictly  by  the  explosion  source  function  in  the 
water  and  should  be  insensitive  to  the  water-bottom  transfer  function  and  the  propaga¬ 
tion-path  effects  of  the  earth.  In  this  study,  we  only  model  the  spectral  and  cepstral  char¬ 
acteristics  of  the  source  function,  not  the  actual  relative  excitation  of  the  regional  phases. 

2.1  UNDERWATER  EXPLOSION  MODEL. 

The  water  column  time-domain  source  function,  W(t),  which  should  be  common  to  all 
phases,  can  be  written  as  follows: 

W(t)  =  X(t)*S(t)*R(t)  (2.1) 

X(t)  is  the  primary  and  bubble-pulse  source  function,  S(t)  is  the  surface-reflection  transfer 
function,  and  R(t)  is  the  water  column  reverberation  transfer  function,  and  *  represents 
the  convolution  function.  We  discuss  each  of  these  terms  below. 
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2.1.1  Bubble  Pulse 


As  mentioned  above,  strong  modulations  observed  in  the  spectra  and  the  sharp  consistent 
peaks  in  the  cepstra  of  these  events  appear  to  result  from  the  interference  of  the  primary 
shock  and  bubble  pulses  in  the  water.  The  concept  of  the  bubble  pulse  is  illustrated  in 
Figure  1. 


Figure  1.  Illustration  of  a  rising  bubble  pulse  from  an  underwater  explosion. 

The  dynamics  of  bubble  pulses,  or  "cavitation,"  have  been  studied  extensively  in  hy¬ 
droacoustics  since  World  War  II.  As  illustrated  in,  after  the  initial  shock  front  is  produced 
by  the  explosion,  there  follow  a  series  of  positive  pressure  pulses  generated  by  the  ex¬ 
panding  and  contracting  gas  sphere  that  rises  to  the  surface.  The  amplitude  of  the  pulses 
decays  steadily  as  the  energy  in  the  expanding  gas  spheres  dissipates  and  as  the  bubble 
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rises  to  the  surface.  Accordingly,  we  model  the  primary  pressure  and  bubble  function 
term  as  follows: 


X(t)  =  Pit )  +  +  y2b(t  -  t2  )  (2.2) 

P(t)  represents  the  primary  pressure  pulse  time  function  in  the  water,  b(t)  is  the  bubble 
pulse  source  function,  and  T2  are  the  first  and  second  bubble  pulse  delay  times.  The 
time  delay  between  the  explosion  pressure  pulse  and  the  first  bubble  pulse  has  been  de¬ 
rived  from  theory  (Willis,  1941)  and  verified  by  a  number  of  studies  (Cole,  1948)  to  be  of 
the  form 


Tl  = 


Kw 


1/3 


(d+ 10) 


5/6 


(2.3) 


where  A'  is  a  proportionality  constant  that  depends  on  the  type  of  explosive,  w  is  the 
charge  size  in  kilograms,  and  d  is  the  depth  of  the  explosions  in  meters.  For  TNT,  K  is 
equal  to  2.1 1.  1  For  the  second  bubble  pulse,  if  present,  we  assumed  Z2=a  Tj,  where  a  is  a 
multiplicative  factor  on  the  order  of  1.8  to  2.0.  These  bubble  pulses  can  continue  until  the 
bubble  rises  and  breaks  the  surface.  For  our  models,  we  allowed  only  one  or  two  bubble 
pulses.  The  multiplicative  amplitude  factor  y in  (2.2)  was  set  to  0.4,  which  was  required 
to  produce  the  spectral  modulations  we  have  observed. 


Figure  2  shows  a  nomogram  of  water  pressure  and  bubble  pulse  frequencies  as  a  function 
of  depth  and  yield.  This  is  the  pressure  nomogram  currently  being  used  by  Air  Force 
Technical  Applications  Center  (AFT AC)  to  estimate  hydroacoustic  pressure  levels  and 
bubble  pulse  frequencies,  and  in  turn,  the  yields  and  depths  of  underwater  explosions  re- 


1  Note  that  in  the  original  paper  of  Baumgardt  and  Der  (1 998),  K  was  erroneously  stated  to  be  equal  to  1.1, 

a  typographical  error.  However,  the  actual  calculations  in  the  paper  were  make  with  the  correct  value. 
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corded  on  a  hydroacoustic  sensor  at  a  distance  of  35°.  For  our  study,  we  are  more  inter¬ 
ested  in  the  bubble-pulse  frequency  since  pressure  levels  are  irrelevant  to  the  analysis  of 
seismic  data.  Bubble-pulse  frequency,  which  is  the  inverse  of  Tj  given  in  (2.1),  is  plotted 
as  dashed  lines.  The  curves  indicate  that  bubble  pulse  frequency  increases  as  a  function  of 
increasing  depth  and  decreasing  yield.  Note  that,  for  continental-shelf  water  depths  of 
about  100  to  1000  ft,  and  yields  between  1  and  100  tons,  the  bubble-pulse  frequencies  are 
in  the  range  of  0.5  to  8  FIz,  well  within  the  short-period  20  Hz  nyquist  seismic  bandwidth 
of  the  data  collection  systems  of  stations  in  the  IMS. 


Yield  (tons) 


Figure  2.  Pressure  and  bubble-pulse  frequency  nomogram.  Solid  lines  are  contours 
of  signal  strength  at  a  distance  of  35°  in  units  of  dB  relative  to  1  mPa. 
Dashed  lines  are  contours  of  the  bubble  pulse  frequency  in  Hz. 

The  pressure  function  can  be  expressed  as  a  step  function  multiplied  by  a  decaying  expo¬ 
nential,  as  follows: 
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P{t )  =  //(/)  exp(-atf), 
H(t)=l  for  /  >  0 
H(0  =  0  for  t  <  0. 


(2.4) 


The  bubble  pulse  term  is  represented  as  a  double-sided  exponential,  as  follows: 


b.(t)  =  exp(-/? t-Ti ), 


(2.5) 


for  the  f  th  bubble  pulse  with  delay  time  Tt.  The  exponential  time  constants  in  (2.4)  and 
(2.5),  are  a=  100  and  (3-  50. 

2,1.2  Water  Column  Bounces 


As  mentioned  above,  an  underwater  explosion  produces  a  primary  pulse  and  one  or  more 
bubble  pulses  resulting  in  a  multiple  pulse  wavelet  that  passes  into  the  water  column. 
This  wavelet  may  then  reflect  off  the  free  water-air  interface  and  pass  back  down  into  the 
water  column,  as  illustrated  schematically  in  Figure  3. 
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Figure  3.  (a)  illustration  of  acoustic  to  seismic  conversion  in  the  water  column,  (b) 
Time  domain  pulses  expected  from  underwater  explosion  and  reflection 
from  the  free  surface,  (c)  Cepstrum  of  direct  and  reflected  pulses.  Cep¬ 
strum  has  negative  peak  at  quefrency  of  the  two-way  reflection  time  in 
water  column. 

For  explosions  near  the  continental  slope,  there  may  only  be  one  reflection  because  the 
laterally  heterogeneous  bathymetry  may  cause  the  bottom  reflections  to  be  scattered  and 
subsequent  reflections  not  observed.  The  spectra  for  near  coastal  underwater  events,  stud¬ 
ied  by  Baumgardt  and  Der  (1998),  seemed  to  be  more  consistent  with  just  a  single  reflec¬ 
tion.  More  examples  will  be  shown  later  in  this  report.  However,  the  model  described  in 
this  section  will  include  other  reflections  as  well. 

In  the  Baumgardt  and  Der  (1994)  model  in  (2.1),  water  column  reflections  were  repre¬ 
sented  as  the  transfer  function  of  the  water  column  reflections  convolved  with  the  first 
surface  reflection.  The  expression  for  the  first  surface  reflection,  S(t),  is 
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(2.6) 


r  J 

s  V  c  ) 


where  d  is  the  depth  of  the  water  column,  c  is  the  speed  of  sound  in  water,  which  is  1480 
m/sec  on  average  for  most  oceanic  regions,  rs  is  the  surface  reflection  coefficient,  and  <Tis 
the  kronecker  delta.  The  transfer  function,  R(t),  for  the  successive  water-column  reflec¬ 
tions  is  given  by 


(2.7) 


where  Mis  the  number  of  reverberations  in  the  water  column  and  rb  is  the  water-bottom 
reflection  coefficient.  We  typically  model  up  to  M—  5  reverberations  in  the  water  column 
after  the  first  surface  reflection. 

Generally,  the  value  of  the  bottom  reflection  coefficient,  rb,  depends  upon  the  water  bot¬ 
tom  geology.  Baumgardt  and  Der  (1998)  found  that  values  of  0.2  and  0.3  fit  most  of  he 
explosions  that  occurred  in  the  Baltic  and  North  Seas.  The  surface  reflection  coefficient, 
rs,  physically  should  be  about  -0.9.  However,  as  we  will  show  later  in  the  report,  we  have 
sometimes  had  to  assumed  lower  values  to  fit  observed  cepstra. 

2.1.3  Synthetic  Cepstra 

Cepstral  analysis  has  proven  to  be  a  very  effective  method  for  characterizing  multiple 
pulse  wavelets,  like  the  ones  produced  by  this  model  of  underwater  explosions.  We  adopt 
the  definition  of  the  “real  cepstrum”,  or  the  rceps  function  from  Matlab  G  Signal  Process¬ 
ing  Toolbox  represented  as  follows: 


RCEPS{x(t)}  =  IFFT{logl0ABS[FFT(x(t)]} 


(2.8) 


Where  FFT  and  IFFT  refer  to  Fast  Fourier  Transform  and  Inverse  Fast  Fourier  transform, 
respectively.  Taking  the  log[0  of  the  spectrum  prior  to  taking  the  next  Fourier  transform 
effectively  whitens  the  spectrum.  Thus,  the  cepstrum  is  just  a  whitened  autocorrelation 
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function  and  the  peaks  that  appear  in  the  cepstrum  are  correlation  peaks. 

In  effect,  the  model  described  above  represents  the  correlation  structure  of  the  waveform 
produced  by  the  source.  As  shown  in  Figure  3,  it  provides  direct  insight  into  the  charac¬ 
teristics  of  the  source  itself  with  no  need  for  more  complex  modeling  of  laterally  hetero¬ 
geneous  earth  structure.  It  is  not  really  necessary  when  modeling  cepstra  to  include  the 
effects  due  to  source  coupling  and  propagation  path  because  only  the  correlation  of  pulses 
is  modeled,  not  the  pulses  themselves.  The  primary  argument  of  this  study  is  that  it  is 
much  easier  and  more  feasible  to  model  cepstra  than  waveforms  or  spectra  since  the  cor¬ 
relation  structure  is  easier  to  model  than  the  waveforms  themselves.  The  proof  of  this  is 
whether  or  not  we  can  adequately  match  synthetic  cepstra  to  observed  cepstra. 

Also,  note  that  the  model  in  expressions  (2.1),  (2.2),  (2.6)  and  (2.7)  will  produce  cepstra 
consisting  of  both  positive  and  negative  peaks,  assuming  that  the  primary  and  secondary 
pulses  from  the  explosion  pressure  waves,  bubble  pulses,  and  water  column  reflections, 
are  correlated.  If  the  two  pulses  are  both  the  same  polarity,  as  they  are  in  the  case  of  the 
primary  pulse  and  the  bubble  pulse,  the  result  in  the  cepstrum  will  be  a  single  positive 
correlation  pulse.  However,  if  there  is  a  free  surface  reflection,  the  primary  and  reflection 
pulse  will  have  opposite  polarity,  as  illustrated  in  Figure  3(b).  Thus,  the  cepstrum,  given 
in,  will  have  a  negative  correlation  peak.  The  composite  cepstrum  will  consist  of  both 
positive  and  negative  peaks,  and  by  matching  the  polarities  of  these  peaks,  we  can  sepa¬ 
rate  the  pulses  produced  by  bubble  pulses  from  those  produced  by  the  surface  reflection. 

2.2  APPLICATION  TO  NOVEMBER  17, 1998  EVENT  NEAR  TROMSO, 
NORWAY. 

Baumgardt  and  Der  (1998)  showed  a  number  of  examples  of  cepstra  with  both  positive 
and  negative  peaks,  presumably  produced  by  bubble  pulses  and  surface  reflections,  re¬ 
spectively.  In  our  Annual  Report  (Baumgardt,  1 999),  we  described  methods  for  matching 
synthetic  cepstra  to  actual  data  cepstra,  demonstrated  on  presumed  explosions  in  the  Gulf 
of  Bothnia  and  in  the  Baltic  Sea.  Here,  we  summarize  the  application  of  the  technique  to 
a  presumed  underwater  event  that  occurred  on  November  17,  1998  near  Tromso,  Norway. 
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2.2.1  Source  and  Waveform  Characteristics 


The  event  occurred  off  the  coast  of  Norway  that  appeared  in  the  Reviewed  Event  Bulletin 
(REB)  for  the  PIDC  with  the  following  event  parameters: 

EVENT  20230910 


Date  Time 

Latitude 

Longitude 

Depth 

Ndef 

Magi 

#  of  Stations 

1998/11/25  11:34:29.5 

69.0645 

16.4273 

0.0  f 

12 

ML  3.5 

4 

This  event  was  not  large  enough  to  be  retained  by  the  REB  at  the  time,  and  so  the  data 
was  not  available  through  the  PIDC.  Data  was  requested  for  the  event  from  NORSAR  for 
the  four  regional  arrays  that  recorded  the  event.  The  map  of  the  propagation  paths  from 
the  event  to  the  arrays  and  a  record  section  of  the  data  are  shown  in  Figure  4. 


Figure  4.  (Left)  Map  showing  propagation  paths  from  the  primary  stations  re¬ 
cording  the  November  25,  1998  event.  The  station  symbols  are  ARCES 
(ARAO),  FINES  (FIC1),  Hagfors  (HFSA1),  and  Kevo  (KEV).  (Right)  Re¬ 
cord  section  of  array  w  aveforms  obtained  from  NORSAR  for  the  event. 

The  regional  phases  recorded  from  the  event  are  indicated  on  the  record  section  on  the 
right.  Regional  phases  Pn,  Sn,  and  Lg  were  observed  across  the  4  stations  ARCES 
(ARAO),  Kevo  (KEV),  FINES  (FIC1),  and  Hagfors  (HFSA1).  The  best  signal-to-noise 
ratios  were  recorded  at  the  ARCES  array  (ARAO)  and  KEV. 

We  relocated  the  event  using  the  travel  times  for  the  regional  phase  picks  made  on  the 
waveforms  shown  in  The  travel  time  models,  1ASPEI91  (Kennett,  1991)  and  the  Baltic 
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Shield  model  of  Bondar  and  Ryaboy  (1997)  were  used  to  relocate  the  event.  The  picks  in 
Figure  4  (right)  are  close  to  the  ones  in  the  REB  except  that  the  NORSAR  picks  were  not 
included.  We  were  not  able  to  see  large  enough  signals  at  NORSAR  to  make  picks  on 
individual  channels  although  the  REB  did  include  picks  for  NORSAR.  It  is  possible  that 
the  REB  picks  were  made  by  the  PIDC  analysis  on  beams,  which  would  enhance  phases 
not  clearly  visible  on  individual  traces. 

Table  1  shows  the  results  of  the  relocation.  The  best  location,  using  the  Baltic  Shield 
model  and  our  phase  picks,  is  shown  in  Figure  5.  IASPEI91  refers  to  the  global  travel 
time  curves  (Kennett,  1991)  that  have  been  used  for  routine  location  in  the  PIDC.  far- 
proved  regional  travel  time  curves  for  this  region  have  been  developed  by  Bondar  and 
Ryaboy  (1997),  which  is  referred  to  as  the  Baltic  Travel  time  curves.  The  term  “local  net¬ 
work”  refers  to  the  phase  picks  in  Figure  4  (b)  and  “IMS  network”  refers  to  the  picks  in 
the  REB. 

Table  1.  Location  Estimates  for  the  November  25, 1998  Presumed  Explosion  Near 
Tromso,  Norway. 


IMS 

IASPEI91 

68.940 

16.583 

2.005 

1.27 

Local 

IASPEI91 

68.935 

16.710 

3679 

0.87 

IMS 

Baltic 

69.052 

16.434 

685 

1.27 

Local 

Baltic 

69.083 

16.465 

1187 

0.96 

This  location  places  the  event  in  amongst  several  groups  of  islands  off  the  coast  of  Nor¬ 
way  and  the  confidence  ellipse  covers  mostly  water-covered  areas.  However,  the  ellipse 
clearly  overlaps  numerous  offshore  islands  although  it  does  not  overlap  the  mainland.  So, 
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refined  location  analysis  in  this  case  cannot  definitively  determine  if  the  event  is  on  land 
or  underwater. 


Figure  5.  Epicenter  and  confidence  ellipse  of  the  relocated  Tromso  event  of  Novem¬ 
ber  28, 1998. 

Figure  6  (top)  shows  a  map  with  the  great  circle  propagation  path  from  the  event  location 
reported  in  the  REB  to  the  ARCES  array  center  element.  Figure  6  (bottom)  shows  part  of 
the  bathymetric/topographic  data  for  the  first  50  km  of  the  path  from  the  source  to  the  re¬ 
ceiver.  The  bathymetry  gives  a  water  depth  of  1 30  m  in  the  vicinity  of  the  source,  al¬ 
though  the  water  depth  varies  from  between  50  to  300  m  along  the  path.  Also,  Figure  6 
shows  that  the  event  was  only  within  about  1 3  to  15  km  of  the  land  of  the  offshore  is- 
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lands,  which  is  less  than  the  REB  epicenter  error  ellipse  long  axis  of  about  21  km.  How¬ 
ever,  as  we  will  show  below,  the  spectral/cepstral  characteristics  of  the  signal  are  consis¬ 
tent  with  a  water  depth  of  1 30  m. 


Distance  from  Source  (Km) 


Figure  6.  (Top)  Great  circle  path  from  the  PIDC  location  of  the  Tromso  event  to  the 
ARCES  array.  (Bottom)  Bathymetric  and  topographic  cross  section  from 
source  to  receiver  showing  water  depth  in  the  region  is  about  130  m. 

2.2.2  Spectral  and  Cepstral  Analysis 


For  the  analysis  of  this  event,  we  chose  to  focus  primarily  on  the  waveforms  from  the 
ARCES  array,  which  was  closest  to  the  event  and  had  the  best  signal-to-noise  ratios  of  all 
the  stations  that  recorded  it. 

Figure  7(a)  shows  a  bandpass  filter  analysis  of  the  ARCES  waveforms,  with  frequency 
bands  from  0.5  to  2.5  Hz  up  to  10  to  16  Hz.  The  signal-to-noise  ratio  is  quite  high 


14 


through  the  entire  spectral  band.  The  Pg  and  Lg  phases  are  notably  large  throughout  the 
band  from  2  to  7  Hz  and  the  first  arrival  Pn  is  emergent.  Pn/Lg  ratios  are  generally  less 
than  1 ,  except  at  frequencies  above  8  Hz,  indicating  significant  shear- wave  excitation. 
Baumgardt  and  Der  (1998)  showed  that,  under  certain  circumstances,  the  Pn/Lg  or  Pg/Lg 
ratio  might  discriminate  earthquakes  and  explosions  for  screening  purposes.  However,  it 
is  doubtful  that  a  Pn/Lg  or  Pg/Lg  ratios  would  confidently  screen  this  event  out  as  an 
earthquake. 


ARAO  Bandpass  Filtered 
Waveforms 
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Figure  7.  (a)  Bandpass  filtered  waveforms  of  the  ARAO  channel  of  the  ARCES  ar¬ 
ray  from  the  November  25, 1998  event  near  Tromso.  (b)  Array  stacked 
spectra  of  the  phases  recorded  at  ARCES. 

Figure  7(b)  shows  the  array-averaged  spectra  for  the  different  phases  and  background 
noise.  The  spectra  are  computed  on  each  individual  ARCES  vertical-component  channel 
by  windowing  the  phases,  starting  at  the  time  picks  shown  in  Figure  7(a),  for  20  seconds 
and  applying  a  Parzen  window.  The  spectra  for  each  channel  are  then  averaged  across  the 
array,  which  produces  much  smoother  spectra  than  the  individual  channel  spectra. 

These  figures  show  spectral  modulations  or  “scalloping”  in  all  phases  with  a  strong  spec¬ 
tral  peak  at  about  3  Hz.  None  of  these  features  appears  in  the  noise.  Spectra  of  this  kind 
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have  been  shown  to  be  caused  by  source  multiplicity,  either  due  to  ripple  fire  in  mine 
blasts  (Baumgardt  and  Ziegler,  1988)  or  by  water  column  reverberations  and  bubble 
pulses  in  underwater  explosions  (Baumgardt  and  Der,  1998).  Since  the  event  occurs  off 
the  coast,  it  is  likely  to  be  an  underwater  explosion,  as  can  be  shown  by  cepstral  analysis. 

The  spectra  in  Figure  7(b)  clearly  have  a  frequency-dependent  trend  due  to  the  combined 
effect  of  the  instrument  response,  the  source  spectral  trend,  and  the  effects  of  the  propaga¬ 
tion  path.  Whitening  the  spectra  would  remove  the  trends  due  to  these  effects  and  en¬ 
hance  the  effect  of  the  spectral  modulations  in  the  cepstum.  We  have  performed  a  number 
of  experiments  with  whitening  the  spectra  for  this  purpose,  hi  our  earlier  research 
(Baumgardt  and  Ziegler,  1988;  Baumgardt  and  Der,  1998),  we  would  first  remove  the  in¬ 
strument  response,  fit  a  polynomial  function  to  the  residual  trend,  and  then  remove  the 
polynomial  trend.  However,  deconvolving  the  instrument  response  sometimes  causes  a 
blowup  in  the  spectrum  at  low  frequencies  that  can  cause  unwanted  processing  artifacts  in 
the  cepstra.  In  our  current  study,  we  have  chosen  to  only  fit  the  polynomial  trend  and  re¬ 
move  it,  skipping  the  step  of  removing  the  instrument  response,  and  have  found  that  this 
procedure  produces  reasonably  stable  cepstra  without  the  low-frequency  blowup  problem. 
This  was  discussed  in  more  detail  in  our  earlier  report  (Baumgardt,  1999). 

The  expression  for  the  cepstrum  can  be  written  as  follows: 

RCEP{x(t)}  =  FFT{  log10  ABSiitf)]}  (2.9) 

where  FFT  stands  for  “fast  fourier  transform”  and  ABSfxJ  refers  to  the  whitened  ampli¬ 
tude  spectrum.  This  expression  is  similar  to  expression  (2.8)  in  the  previous  section  for 
the  theoretical  cepstra  except  that  the  spectra  have  been  whitened  by  removing  the  poly¬ 
nomial  trend  prior  to  taking  the  second  FFT.  In  the  case  of  theoretical  cepstra,  no  whiten¬ 
ing  is  applied  to  the  spectrum,  other  than  the  logio,  before  taking  the  second  FFT. 

For  the  regional  arrays,  cepstra  can  be  computed  for  each  of  the  phases  and  then  averaged 
across  the  array  channels.  This  averaging  accentuates  the  peaks  and  troughs  in  the  cepstra 
caused  by  real  source  effects,  since  they  show  up  on  each  channel,  and  diminish  those  due 
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to  site  effects.  Comparisons  of  single  channel  and  array  averaged  cepstra  have  shown  that 
the  array  averaged  cepstra  have  much  sharper  peaks  and  troughs  than  the  single  channel 
cepstra.  Moreover,  further  enhancement  of  the  peaks  and  troughs  can  be  attained  by  aver¬ 
aging  or  stacking  the  array-averaged  cepstra  for  each  phase  to  form  a  composite,  stack 
cepstrum. 

Figure  8  shows  cepstra  for  the  three  phases,  the  stack  of  the  three  cepstra,  and  the  cepstra 
of  the  background  noise.  The  resultant  cepstra,  using  (2.9),  are  a  function  of  delay  time 
between  the  multiple  pulses,  called  “quefrency”  in  units  of  time  (sec),  and  are  signed.  A 
cepstrum  for  the  noise  spectrum  is  also  computed  in  order  to  check  for  cepstral  features 
that  may  be  due  to  processing  artifact  rather  than  source  effects.  We  look  for  cepstral 
peaks  in  the  individual  phase  cepstra  that  do  not  appear  in  the  noise  cepstra.  The  cepstra 
in  have  been  shifted  for  display  purposes  and  the  absolute  amplitudes  are  arbitrary.  Posi¬ 
tive  cepstral  peaks  indicate  pulses  that  are  correlated  with  the  same  polarity  and  negative 
peaks,  or  troughs,  are  produced  by  pulses  with  reverse  polarity.  As  shown  by  Baumgardt 
and  Der  (1998),  positive  peaks  are  caused  by  bubble  pulses  with  positive  polarity  and 
negative  peaks  are  produced  by  reflections  of  acoustic  waves  off  the  water-air  interface. 
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Figure  8.  Shifted  cepstra  computed  for  the  spectra  shown  in  Figure  7(b). 

2.2.3  Synthetic  Cepstra 

The  modeling  techniques,  described  by  Baumgardt  and  Der  (1998)  and  in  the  previous 
section,  are  illustrated  in  Figure  9.  We  compute  wavelets  for  the  bubble  pulse,  shown  on 
top,  and  water  column  reverberations,  shown  as  the  second  trace,  for  a  given  assumed 
yield  and  depth  in  the  water.  Other  parameters  that  must  be  specified  include  average  wa¬ 
ter  column  depth,  that  can  be  obtained  from  the  bathymetry,  and  the  surface  and  water 
bottom  reflection  coefficients,  which  we  assume  to  be  -0.9  and  0.3,  respectively.  As  in 
the  case  of  water  depth,  these  parameters  can  be  estimated  with  some  degree  of  accuracy, 
but  it  is  also  conceivable  that  they  can  be  inferred  from  the  data.  We  assume  that  the 
acoustic  waves  in  the  water  reflect  from  first-order  discontinuities  at  the  bottom  and  the 
surface. 
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Figure  9.  Method  for  computing  the  cepstral  model  for  matching  observed  cepstra. 

The  final  synthetic  cepstmm  is  the  convolution  of  the  source  plus  bubble  pulse  wavelet 
and  reflection  wavelet.  We  then  compute  the  cepstrum  of  this  composite  wavelet,  which 
is  shown  at  the  bottom  of  Figure  9.  We  also  show  a  minimum  phase  wavelet  reconstruc¬ 
tion  from  the  cepstrum,  which  is  identical  to  the  original  wavelet  since  the  wavelets  were 
constructed  minimum  phase.2 

2.2.4  Cepstral  Inversion 

Finally,  we  seek  a  synthetic  cepstrum,  like  that  at  the  bottom  of  Figure  9  that  matches  the 
observed  cepstrum.  For  this  purpose,  we  compute  the  correlation  coefficient  between  the 


2  The  assumption  of  "minimum  phase”  means  that  the  amplitudes  of  the  pulses  in  the  wavelet  decrease  with 
time,  a  reasonable  assumption  for  bubble  pulses  and  surface  reflections.  However,  there  is  no  way  from  the 
cepstrums  alone  to  determine  if  the  signals  are  minimum  phase  or  not. 
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stacked  observed  cepstrum  and  the  synthetic  cepstrum  and  seek  a  cepstrum  that  provides 
the  maximum  correlation. 

We  have  first  considered  a  simple  exhaustive  search  method,  where  the  parameter  space 
of  blast  yield  and  depth  in  the  water  is  gridded  and  synthetic  cepstra  are  computed  for 
each  grid  point.  The  method  correlates  the  synthetic  cepstrum  at  each  grid  point  with  the 
observed  cepstrum,  compute  a  correlation  coefficient  surface  over  the  depth/yield  pa¬ 
rameter  space,  and  look  for  the  peaks  in  the  depth/yield  surface. 

Figure  10  shows  the  resultant  correlation  coefficient  surface  plotted  as  a  3D  surface  (left) 
and  as  a  contour  plot  (right).  Both  displays  have  a  distinct  peak  in  the  correlation,  with  a 
maximum  correlation  coefficient  of  0.73,  corresponding  to  a  water  depth  of  1 30  m  and 
yield  of  1800  kg  or  about  1.98  tons.  It  is  evident  in  that  depth  is  well  constrained  but  yield 
is  less  well  constrained.  The  yield  is  controlled  exclusively  by  the  bubble  pulse  delay 
time,  although  the  bubble  pulse  delay  is  also  a  function  of  depth.  However,  depth  is  con¬ 
strained  also  by  the  negative  peak  delay  time,  corresponding  to  the  two-way  water- 
column  reflection  time. 


Figure  10.  Synthetic  cepstra  computed  for  a  range  of  yield  (1000  to  3000  kg)  and 

depth  (0  to  200  m)  and  each  correlated  with  the  observed  cepstra.  Surface 
plot  (left)  and  contour  plot  (right)  showing  the  maximum  peak  that  gives 
the  best  estimate  of  the  depth  and  yield  of  the  explosion. 


20 


The  search  algorithm  found  a  depth  of  130  m  based  on  the  0.2-second  surface  reflection 
time.  Note  that  this  would  be  the  depth  of  the  explosion  below  the  water  surface.  It  is  no¬ 
table  that  the  inferred  explosion  depth  of  130  m  is  close  to  the  actual  water  depth  of  the 
source  location  as  inferred  from  bathymetry,  shown  in  the  cross  section  in  Figure  6. 

This  suggests  that  the  explosion  may  have  occurred  near  the  bottom  or  on  the  bottom. 
However,  that  knowledge  of  the  water  depth  is  not  actually  required  to  infer  the  explosion 
depth  in  the  water,  since  the  reflection  from  the  surface  is  by  far  the  strongest  signal.  We 
actually  assumed  a  water  depth  of  225  m  in  our  model.  Reflections  from  the  bottom  may, 
in  principle,  be  observed  and  reveal  water  depth,  but  they  are  very  weak,  as  is  evident  in 
Figure  9. 

Generally,  we  have  not  observed  strong  bottom  reflections  in  the  observed  cepstra,  which 
indicates  that  the  bottom  reflection  coefficient  is  very  small.  Small-reflection  coefficients 
may  be  caused  by  roughness  of  the  bottom  itself.  Moreover,  transitional  ocean  bottom 
layers  may  cause  frequency-dependent  reflection  coefficients  that  distort  the  bottom  re¬ 
flections  such  that  they  do  not  correlate  with  the  primary  pulses  and  produce  strong  cep- 
stral  peaks.  However,  water  depths  at  given  locations  in  the  ocean  are  usually  well  known 
and  need  not  be  inferred  by  this  method.  However,  knowing  the  maximum  water  depth  is 
useful  because  it  constrains  the  required  search  space. 

The  sensitivity  to  yield  is  shown  in  Figure  1 1,  where  the  correlation  coefficient  is  plotted 
as  a  function  of  yield  for  the  depth  of  1 30  m.  The  maximum  peak  was  at  1 800  kg,  but  an¬ 
other  peak  of  comparable  strength  appears  at  closer  to  2000  kg.  However,  the  overall 
broad  peak  is  about  200  or  300  kg  in  width,  which  provides  an  estimate  of  the  resolution 
of  the  method  to  yield.  Thus,  a  better  estimate  of  yield  might  be  1900  kg  ±  200  kg. 
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Figure  11.  Sensitivity  of  the  correlation  method  for  estimating  underwater  explosion 
yield  by  matching  synthetic  and  observed  cepstra. 

Finally,  Figure  12  shows  a  direct  comparison  of  the  observed  stack  cepstrum  and  the  best¬ 
matching  synthetic  cepstrum.  This  comparison  shows  that,  although  the  match  is  not  per¬ 
fect,  reasonably  good  matches  of  the  essential  features,  the  negative  peaks  of  the  water 
column  reflection  and  the  bubble  pulses,  have  been  attained.  The  primary  features  that 
drive  the  match  correlation  is  the  position  of  the  negative  surface-reflection  peak  and  the 
bubble  pulse  peaks,  three  of  which  can  be  seen  in  this  case. 

We  can  also  check  the  yield  with  the  AFTAC  pressure  nomogram  in  Figure  2.  The  bubble 
pulse  period  for  the  best  fitting  cepstrum  in  Figure  12  is  about  0.42  sec,  or  a  bubble-pulse 
frequency  of  2.4  Hz.  (In  spectral  terms,  the  spectral  scalloping  that  gives  rise  to  the  cep- 
stral  peak  has  modulation  period  of  2.4  Hz.)  Given  a  depth  of  130  m,  or  426  ft,  Figure  1 
gives  a  yield  of  a  just  above  1  ton,  that  is  consistent  with  the  1 .98  ton  value  obtained  by 
cepstral  fitting  search  algorithm. 
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Figure  12.  Match  of  the  synthetic  cepstra  (dotted)  to  the  observed  stacked  cepstra 
the  best  fitting  underwater  explosion  model. 

Thus,  this  analysis  strongly  suggests  that  this  event  must  be  an  underwater  explosion  of 
close  to  2  tons  to  explain  the  spectra  and  cepstra.  The  event  was  not  identified  by  the 
PIDC  as  an  underwater  explosion,  and  in  fact,  the  IMS  hydroacoustic  network  did  not  de¬ 
tect  the  event  because  propagation  paths  were  blocked  to  all  hydroacoustic  stations. 

It  might  be  argued  that  this  is  only  a  small  conventional  explosion  of  no  importance  and  if 
it  had  been  a  nuclear  blast  of  much  larger  yield,  the  IMS  hydroacoustic  network  would 
have  easily  detected  it.  From  Figure  2,  the  pressure  level  at  130  m  for  a  1-ton  shot  would 
have  been  between  70  and  75  mPa  at  a  distance  of  35°.  If  the  explosion  had  been  10,000 
tons,  or  10  Kt,  consistent  with  a  small  nuclear  blast,  the  pressure  levels  at  35°  from  Figure 
1  would  only  be  between  85  and  90  mPa.  This  increase  in  pressure  level  for  a  two  order- 
of-magnitude  increase  of  yield,  may  not  be  able  to  overcome  the  blockage  effects  in  order 
to  be  detected  at  an  IMS  hydroacoustic  station.  Thus,  this  kind  of  analysis  should  also  be 
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applied  seismic  data  as  well  as  hydroacoustic  data  for  offshore  events  near  the  coast. 

Several  other  examples  of  presumed  underwater  chemical  blasts  were  described  in  our 
earlier  report  (Baumgardt,  1 999).  However,  we  had  no  actual  ground  truth  on  the  events, 
although  the  data  was  compelling  that  they  were  actually  explosions.  In  order  to  validate 
that  the  technique,  we  need  to  test  the  method  on  actual  underwater  explosions  for  which 
ground-truth  information  is  available.  In  the  next  section,  we  discuss  the  application  of 
the  method  to  the  DTRA  sponsored  calibration  explosions  detonated  in  the  Dead  Sea  in 
November  of  1999. 
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SECTION  3 

CEPSTRAL  CHARACTERZATION  OF  THE  NOVEMBER  8, 10,  and  11 
CALIBRATED  UNDERWATER  EXPLOSIONS  IN  THE  DEAD  SEA 


In  November  of  1999,  under  the  sponsorship  of  DTRA,  three  calibration  chemical  explo¬ 
sions  were  detonated  in  the  Dead  Sea  by  the  Geophysical  Institute  of  Israel.  The  main  ob¬ 
jective  of  the  detonations  was  to  calibrate  seismic  travel  times  at  local  and  regional  dis¬ 
tances  in  the  region.  Each  of  the  three  explosions  are  at  three  different  yields  and  at  the 
same  depth  of  70  m.  The  event  parameters  for  the  explosions  were  (Gitterman  et  al,  1999; 
Gittermann  and  Shapira,  2000)  are  given  in  Table  2. 

Table  2.  Event  Parameters  for  the  three  Dead  Sea  Calibration  Events. 


1999/11/08 

13:00:00.33 

31.5330 

35.4406 

70  m 

3.1 

- 

500 

1999/11/10 

13:59:58.21 

31.5338 

35.4400 

70  m 

3.6 

- 

2060 

1999/11/11 

15:00:00.79 

31.5336 

35.4413 

70  m 

3.9 

4.1 

5000 

The  detonations  were  timed  with  GPS  with  an  accuracy  of  5  milliseconds.  The  location 
accuracy  was  on  the  order  of  50  m. 


Our  analyses  of  these  events  allowed  the  testing  of  our  developed  algorithm  to  known, 
calibrated  explosions,  hi  this  section,  we  describe  the  cepstral  analysis  and  inversion  of 
these  events  in  order  to  validate  the  accuracy  of  the  method  for  estimating  explosion 
yield. 

For  our  analysis,  we  used  data  recorded  at  regional  IMS  stations  available  to  us  from  the 
PIDC.  The  locations  of  the  regional  Israeli  IMS  stations  at  Eilat  (EIL)  and  Meron  (MRNI) 
that  recorded  the  Dead  Sea  explosions  are  shown  in  Figure  13.  The  figure  shows  a  map 
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with  the  great  circle  propagation  paths  from  the  locations  of  the  explosions  to  the  MRNI 
and  EIL  sites. 


Figure  13.  Great  Circle  Paths  from  Dead  Sea  Explosions  to  MRNI  and  EIL. 


3.1  DEAD  SEA  EXPLOSION  WAVEFORMS. 

Figure  14  shows  three-component  MRNI  waveforms  for  the  November  9,  1999  explo¬ 
sion.  Data  from  the  EIL  station  was  not  available  for  this  event.  Waveforms  recorded  on 
the  EIL  and  MRNI  stations,  for  the  November  10  and  1 1,  1999  explosions  are  shown  in 
Figure  15  and  Figure  16,  respectively.  Figure  17  shows  the  waveforms  for  both  the  EIL 
and  MRNI  stations  for  the  November  11,  1999  explosion.  The  regional  phase  picks  made 
by  the  PICD  analysts  are  shown  on  each  of  the  waveforms. 

Figure  17  compares  the  Eilat,  Israel  (EIL)  recordings  of  November  19  and  November  11, 
1999  Dead  Sea  explosions.  Figure  18  compares  the  Meron,  Israel  (MRNI)  recordings  of 
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November  8,  10,  and  1 1,  1999.  These  comparisons  show  that  the  relative  amplitudes  of 
the  waveforms  approximately  correlate  with  the  differing  explosion  yields. 


November  8,  1999  Dead  Sea  Event  -  MRNI 


Figure  14.  November  8,  1999  Explosion:  Eilat,  Israel  (EIL),  Three-Component 
Phase  Picks. 
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Figure  15.  November  10, 1999  Explosion:  (left)  Eilat,  Israel  (EIL)  Three-Com¬ 
ponent  Phase  Picks;  (right)  Meron,  Israel  (MRNI)  Three-Component 
Phase  Picks. 


Figure  16.  November  11, 1999  Explosion:  (left)  Eilat,  Israel  (EIL)  Three-Component 
Phase  Picks;  (right)  Meron,  Israel  (MRNI)  Three-Component  Phase 
Picks. 
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Comparison  of  EIL  Vertical  Recordings  for  Two  Evenls  -  PIDC  Phase  Picks 


2000  Kg  Explosion 


Figure  17.  Comparison  of  EIL  Vertical  Recordings  for  November  10,  and  November 
11,  1999  Explosions. 


Comparison  of  MRNI  Vertical  Recordings  for  Three  Events  -  PIDC  Phase  Picks 
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Figure  18.  Comparison  of  MRNI  Vertical  Recordings  of  November  8, 10,  and  11 
Explosions. 
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3.2  SPECTRAL/CEPSTRAL  ANALYSIS. 


We  use  the  same  method  for  computing  spectra  and  cepstra  that  was  described  in  the  ear¬ 
lier  discussion  of  the  Tromso.  Figure  19  shows  the  vertical-component  spectra  for  the  dif¬ 
ferent  phases  and  background  noise,  for  the  phases  recorded  at  the  for  the  three  Dead  Sea 
underwater  explosions.  The  spectra  were  computed  on  each  of  the  individual  three- 
component  channels  by  windowing  the  phases,  starting  at  the  time  picks  shown  in  Figure 
17  and  Figure  18,  using  a  20-second  Parzen  window.  Since  no  array  data  was  available 
for  these  events,  the  spectra  for  each  channel  were  averaged  across  the  three  component 
waveforms,  which  produce  much  smoother  spectra  than  the  individual  channel  spectra. 
The  spectra  have  very  distinct  spectral  modulations  or  “scalloping”  in  all  phases  with  a 
strong  spectral  peak  at  about  2.5  Hz.  None  of  these  features  appears  in  the  noise. 


Figure  19.  Averaged  spectra  for  each  of  the  recordings  of  the  Dead  Sea  underwater 
blasts. 
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Figure  19  also  shows  the  “usable  bandwidth”  of  each  of  the  signal.  The  “usable  band¬ 
width”  is  the  portion  of  the  spectrum  where  the  phase  spectral  ampl  itudes  are  above  the 
noise  levels.  Within  the  usable  bandwidth  in  each  case,  the  spectral  modulations  are  evi¬ 
dent,  although  they  are  stronger  for  the  larger  explosions  on  November  10  and  1 1  than  for 
the  November  8  explosion.  It  is  obvious  from  Figure  19  that  the  “usable  bandwidth”  is 
considerably  smaller  than  the  overall  bandwidth  of  the  data.  In  order  to  apply  the  same 
cepstral  analysis  as  that  described  in  the  previous  section  for  synthetic  data,  spectral  shap¬ 
ing  is  required. 


We  use  the  same  technique  for  computing  cepstra  of  the  Dead  Sea  events  as  we  used  for 
the  Tromso  events,  which  is  expression  (2.9).  Figure  20  shows  the  spectra  (left)  and  cep¬ 
stra  (right)  for  the  MRNI  recordings  of  the  November  8,  1999  explosion.  In  this  figure, 
both  the  spectra  and  the  cepstra  have  been  shifted  for  display  purposes.  This  explosion 
was  the  smallest  of  the  three  and  has  the  weakest  spectral  modulations,  and  we  only  have 
data  from  this  one  station.  Figure  21  and  Figure  22  show  the  composite  spectra  for  the 
November  10  and  1 1  explosions,  respectively.  For  these  two  events,  the  spectra  for  all 
channels  of  MRNI  and  EIL  have  been  averaged. 
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Figure  20.  Spectra  and  cepstra  for  the  November  8, 1999  explosion. 
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Figure  21.  Averaged  MRNI  and  EIL  spectra  and  cepstra  for  the  November  10,  1999 
explosion. 


EIL  and  MRNI  Shifted  Spectra 


EIL  and  MRNI  Shifted  Cepstra 


-  Ns 

-  Pn 

PQ 

Lg 

-  Stack 


Figure  22.  Averaged  MRNI  and  EIL  spectra  and  cepstra  for  the  November  11,  1999 
explosion. 

These  plots  show  the  spectra  and  cepstra  for  each  phase  that  was  identified  and  the  back¬ 
ground  noise  sample  before  the  Pn,  and  a  composite  stack  of  all  the  phase  cepstra  for 
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each  event.  The  strongest  modulations  were  observed  for  the  November  10  and  1 1  events, 
as  might  be  expected  since  they  are  the  largest.  The  larger  the  explosion,  the  higher  the 
signal  relative  to  the  noise  in  the  usable  bandwidth,  the  stronger  the  spectral  modulations, 
and  the  larger  the  cepstral  peaks.  Also,  the  period  of  the  primary  peak,  which  as  associate 
with  the  bubble  pulse,  scales  with  the  size  of  the  explosion.  The  November  10  event  ap¬ 
pears  to  have  peaks  that  may  either  be  caused  by  two  bubble  pulses  or  multiple  explo¬ 
sions.  This  event  did  have  a  small  secondary  charge  but  it  was  only  about  60  kg.  We  now 
apply  the  cepstral  model  and  inversion  algorithm  on  the  composite  stack  cepstrum  of  each 
event. 

3.3  CEPSTRAL  MODELING  AND  INVERSION. 

We  have  first  considered  a  simple  exhaustive  search  method,  where  the  parameter  space 
of  blast  yield  and  depth  in  the  water  is  gridded  and  synthetic  cepstra  are  computed  for 
each  grid  point,  as  we  did  for  the  Tromso  event  in  the  previous  section.  However,  for  the 
match  statistic,  we  have  used  the  LI  norm  rather  than  the  correlation  coefficient  to  find 
the  best  matching  synthetic  cepstrum  to  the  observed  cepstrum.  The  LI  norm  is  the  abso¬ 
lute  value  of  the  difference  between  the  observed  and  synthetic  cepstrum,  which  would  be 
zero  for  a  perfect  match.  Thus,  in  the  exhaustive  search,  we  seek  the  parameters  of  explo¬ 
sion  depth  and  yield  that  maximize  the  1.0  -  LI  norm. 

The  result  of  the  exhaustive  search  solution  to  the  fit  to  the  November  8  explosion  is 
shown  in  Figure  23.  Although  a  peak  is  found  that  gives  a  yield  of  650  kg,  close  to  the 
expected  value  of  500  kg,  the  peak  is  not  well  resolved. 
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Contour  Plot  of  Depth  vs  Yield  3D  Plot  of  Depth,  Yield,  and  1.0-  II norm 


Inversion  Yield  =  650kg,  Depth  =  80m  Actual  Yield  =  500kg,  Depth  =  70m 


Figure  23.  LI  norm  exhaustive  search  solution  to  the  November  8, 1999  event.  The 

peaks  indicate  the  best  matching  yield  and  depth  on  the  contour  plot  (left) 
and  surface  plot  (right)  for  the  1  -  LI  norm. 

Figure  24  (left)  shows  the  sensitivity  of  the  data  to  this  solution,  where  the  1  -  LI  norm 
coefficient  is  plotted  as  a  function  of  yield  for  a  depth  of  80  m.  The  maximum  peak  was 
650  kg.  The  overall  broad  peak  is  about  140  kg  in  width,  which  provides  an  estimate  of 
the  resolution  of  the  method  to  yield;  thus,  a  better  estimate  of  yield  is  650  kg  ±  140  kg. 
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Figure  24.  November  8,  1999:  (left)  cross  section  of  cepstral  fit  surface;  (right)  com¬ 
parison  of  synthetic  cepstrum  with  stack  MRNI  data  cepstrum. 
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Figure  24  (right)  shows  a  direct  comparison  of  the  observed  stack  cepstrum  and  the  best¬ 
matching  synthetic  cepstrum.  This  comparison  shows  that,  although  the  match  is  not  per¬ 
fect,  the  theoretical  cepstrum  matches  the  overall  shape  of  the  observed  cepstrum,  but  not 
the  details.  The  primary  feature  that  drives  the  match  correlation  is  the  position  of  the 
bubble  pulse  peak.  It  should  be  noted,  however,  that  the  depth  is  not  well  constrained. 
Also,  there  appears  to  be  significant  low-queffency  noise  due  the  processing  artifacts. 


The  solution  for  the  November  10,  1999  event  is  shown  in  Figure  25.  This  solution  pro¬ 
duced  multiple  peaks  in  the  match  surface.  The  maximum  peak  gave  a  yield  of  4000  kg, 
which  is  considerably  larger  than  the  true  yield  of  2000  kg.  However,  it  should  be  noted 
that  the  depth  of  the  estimate  of  the  maximum  peak  is  110  m,  which  is  considerably 
greater  than  the  true  depth  of  70  m.  The  secondary  peaks  correspond  to  lower  depths,  in¬ 
cluding  one  at  about  80m. 


o> 

JC 

"O  ooo 

V 


Contour  Plot  of  Depth  vs  Yield 


Maximum  Peak 


1.0-11  norm  =  .96 


Secondary 

Peaks 


3D  Plot  of  Depth,  Yield,  and  1-llnorm 

Secondary 

o^laximum  Peak 
jjQ-JI  norm  =  .96 


Peaks 


B5  70  75 


»5  tOO  106  110 


Yield  (kg) 


Depth  (m) 

Yield  =  4000  kg,  Depth  =  1 10  m  Actual  Yield  =  2000  kg,  Depth  =  70  m 


Depth  (m) 


Figure  25.  LI  norm  exhaustive  search  solution  to  the  November  10, 1999  event.  The 
peaks  indicate  the  best  matching  yield  and  depth  on  the  contour  plot  (left) 
and  surface  plot  (right)  for  the  1  -  LI  norm. 

Figure  26  shows  the  sensitivity  to  yield  where  the  1  -  LI  Norm  coefficient  is  plotted  as  a 
function  of  yield  for  the  depth  of  80  m.  The  maximum  peak  was  at  1950  kg,  which  is 
closer  to  the  true  value.  The  overall  broad  peak  is  about  700  kg  in  width,  which  gives  a 
yield  estimate,  with  uncertainty,  of  1950  kg  ±  350  kg. 
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Figure  26.  November  10, 1999:  (left)  cross  section  of  cepstral  fit  surface;  (right) 
comparison  of  synthetic  cepstrum  with  stack  EIL  and  MRNI  data  cep- 
strum. 

Figure  26  (right)  shows  a  direct  comparison  of  the  observed  stack  cepstrum  and  the  best¬ 
matching  synthetic  cepstrum  at  depth  of  80  m.  This  comparison  shows  that  a  double¬ 
bubble  pulse  gives  a  good  match  of  the  essential  features.  Thus,  the  complexity  in  the  fit¬ 
ting  surface,  caused  by  the  multiple  peaks,  can  be  explained  by  a  double  bubble-pulse  fea¬ 
ture  in  the  cepstrum.  Again,  the  depth  appears  to  be  poorly  constrained,  with  significant 
low-quefrency  noise  evidently  caused  by  processing  artifact. 


Figure  27  shows  the  solution  for  the  November  11,  1999  explosion  plotted  as  a  3D  con¬ 
tour  (left)  and  surface  plot  (right).  Both  displays  have  a  distinct  peak  in  the  correlation, 
with  a  maximum  1  -  LI  Norm  coefficient  of  0.98,  corresponding  to  a  water  depth  of  80 
m  and  yield  of  4200  kg.  The  secondary  peak  of  .976  is  also  shown,  which  gives  a  yield 
closer  to  the  expected  value  of  5000  kg  and  but  the  depth  is  90  m. 
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Figure  27.  November  11, 1999:  Example  of  exhaustive  search  for  best  matching  syn¬ 
thetic  cepstrum  to  observed  cepstrum  1  -  LI  Norm,  (left)  contour  plot 
(right)  surface  plot. 

The  sensitivity  to  yield  is  shown  in  Figure  28  (left),  where  the  maximum  1  -  L2  norm 
correlation  is  plotted  as  a  function  of  yield  for  the  depth  of  80  m.  The  maximum  peak  was 
at  4200  kg.  However,  the  overall  broad  peak  is  about  800  kg  in  width,  which  gives  an  es¬ 
timate  of  yield  and  uncertainty  of  4200  kg  ±  500  kg.  Figure  28  (right)  shows  a  direct 
comparison  of  the  observed  stack  cepstrum  and  the  best-matching  synthetic  cepstrum. 
Note  that  there  is  a  very  strong  negative  peak,  which  evidently  comes  from  the  surface 
reflection.  This  negative  peak  provides  the  constraint  on  depth. 
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Figure  28.  November  11, 1999:  (left)  cross  section  of  cepstral  fit  surface;  (right) 
comparison  of  synthetic  cepstrum  with  stack  EIL  and  MRN1  data  cep¬ 
strum. 

The  negative  surface  reflection  is  about  0.15  seconds.  We  assumed  a  water  velocity  of 
1771  m/sec,  which  is  higher  than  the  oceanic  value  of  about  1500  m/sec  because  of  the 
high  salinity  of  the  Dead  Sea  (Gitterman  et  al,  1999).  This  time  delay  is  consistent  with  a 
two-way  travel  time  through  a  water  column  of  about  80  m,  which  was  found  by  the 
search  algorithm.  Note  that  this  would  be  the  depth  of  the  explosion  below  the  water  sur¬ 
face.  The  actual  charge  depth  was  70  m.  As  noted  earlier,  knowledge  of  the  water  depth  is 
not  actually  required  to  infer  the  explosion  depth  in  the  water  since  the  reflection  from  the 
surface  is  by  far  the  strongest  signal. 

For  all  the  cepstral  matches  described  above,  we  have  had  to  assume  a  very  low  surface 
reflection  coefficient  of -0.5,  considerable  lower  than  the  expected  value  of -0.90.  This 
lower  value  is  required  to  match  the  negative  trough  amplitude.  A  coefficient  of  -0.9 
would  make  the  synthetic  trough  almost  twice  as  large  as  the  observed  trough.  The  low 
amplitude  of  the  trough  in  the  model  cepstrum  is  evidently  caused  by  low  correlation  of 
the  primary  explosion  pulse  and  the  reflection.  Because  of  the  shallowness  of  the  explo¬ 
sion,  we  are  barely  able  to  resolve  the  reflection  with  the  limited  bandwidth  of  the  data. 
Thus,  this  poor  resolution  may  cause  the  trough  to  be  weaker  than  observed  theoretically. 
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Regardless  of  the  value  assumed  for  the  surface  reflection  coefficient,  only  the  time  of  the 
surface-reflection  trough  in  the  data  needs  to  be  matched  to  the  model,  not  the  amplitude 
itself,  in  order  to  estimate  depth. 

3.4  DOWNHILL  SIMPLEX  SEARCH  ALGORITHM. 

The  exhaustive  grid-search  cepstral  matching  algorithm  that  was  applied  previously  has 
the  disadvantage  of  being  very  time  consuming,  and  in  its  current  implementation,  we 
could  only  vary  the  yield  and  depth  parameter  of  the  model.  A  number  of  other  important 
parameters,  such  as  water  depth,  surface  reflection  coefficient,  and  bottom  reflection  co¬ 
efficients,  had  to  be  fixed.  A  faster  algorithm  that  can  perturb  more  parameters  would  be 
more  desirable. 

Here,  we  have  applied  the  downhill  simplex  method,  originally  developed  by  Nelder  and 
Mead  (1965),  for  multidimensional  search  of  a  parameter  space  to  minimize  a  functional. 
Given  an  A  dimensional  parameter  space,  the  algorithm  searches  the  complex  topography 
of  the  function  surface  in  the  downhill  direction  to  minimize  a  goodness  of  fit  parameter. 
In  our  case,  the  function  is  defined  as  the  match  statistic  (e.g.,  correlation  coefficient  or  1- 
L1  norm),  R[  C0,  Cth(x)],  between  the  observed  cepstrum,  C0,  and  the  theoretical  cep- 
strum,  Cth,  and  x  is  an  A  dimensional  vector  of  parameters  for  the  theoretical  cepstrum. 
We  desire  parameters  for  a  model  that  minimize  the  function  £  =  1  -  R[  C0,  Cth(x)]. 

Very  briefly,  the  downhill  simplex  algorithm  finds  the  minimum  of  the  functional  of  an  A 
dimensional  surface  by  first  defining  a  simplex  with  A+7  points.  Given  a  starting  point 
P„,  the  A  other  points  are  taken  as 


Pi  —  P0  +  A,  e* 

where  the  ej-s  are  A  unit  vectors,  and  where  X  is  a  constant  that  defines  the  guess  values 
initial  length  scale.  The  downhill  simplex  method  then  takes  a  series  of  steps,  moving 
from  the  point  of  the  simplex  where  the  function,  £,  is  largest  through  the  opposite  face  of 
the  simplex  to  a  lower  point.  The  method  tries  to  find  a  valley  and  shrinks  the  simplex  to 
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better  define  the  minimum  in  the  valley. 


In  our  initial  application  of  the  algorithm,  we  allowed  three  parameters  to  be  perturbed  in 
the  model  —  depth,  yield,  and  surface  reflection  coefficient.  The  reason  surface  reflection 
coefficient  was  included  as  a  parameter  is  that  lower  values  than  -0.9  were  needed  to  fit 
the  amplitude  of  the  negative  peak  in  the  cepstrum.  However,  given  that  we  are  using  the 
cepstral  correlation  as  the  goodness-of-fit  functional,  it  is  not  really  necessary  to  have  an 
exact  match  in  the  amplitudes  of  the  peaks  in  order  to  get  a  high  correlation.  Still,  as  an 
experiment  with  the  method  and  to  find  out  what  is  necessary  to  get  the  synthetic  cepstra 
to  look  more  like  the  observed  cepstra,  we  have  allowed  the  surface  reflection  coefficient 
to  vary. 

The  downhill  simplex  algorithm,  originally  developed  by  Nelder  and  Mead  (1965),  com¬ 
putes  a  multidimensional  search  of  depth,  yield,  and  surface  reflection  coefficient  pa¬ 
rameter  space  to  maximize  the  correlation  between  the  synthetic  cepstrum  and  the  data 
cepstrum.  Figure  29  (left)  shows  a  plot  of  the  parameter  space,  consisting  of  yield  and 
depth  on  the  horizontal  axes,  and  surface  reflection  coefficient  on  the  vertical  axis.  Each 
‘x’  symbol  marks  a  point  in  the  parameter  space  that  was  checked  by  the  downhill  sim¬ 
plex  method.  The  arrow  indicates  the  point  that  the  algorithm  found  as  the  minimum  of 
the  LI  norm  match  statistic.  This  is  the  point  where  the  simplex  algorithm  collapses  in  on 
a  minimum  value.  The  right  plot  of  Figure  29  shows  the  best  match,  found  by  the  down¬ 
hill  simplex  optimization  method,  between  the  data  cepstrum  and  the  synthetic  cepstrum. 
The  match  is  good  overall,  although  the  height  of  the  synthetic  bubble  pulse  is  not  as  big 
as  that  of  the  data  cepstrum.  The  downhill  simplex  algorithm  chose  a  yield  of  4200  kg, 
and  a  depth  of  80  m.  The  actual  yield  was  5000  kg,  at  a  depth  of  70  m. 
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Figure  29.  November  11,  1999:  (left)  downhill  simplex  optimization  method;  (right) 
comparison  of  synthetic  cepstrum  with  highest  correlation. 


3.5  CONCLUSIONS. 

The  limited  bandwidth  of  the  IMS  regional  stations  was  not  sufficient  to  resolve  the  70  m 
depth  of  the  explosions.  However,  we  have  shown  that  the  cepstra  for  the  explosions  are 
consistent  with  the  known  depth  of  70  m.  The  November  10,  1999  event  appeared  to  have 
two  bubble  pulses,  or  perhaps  to  detonations.  Lack  of  constraint  on  depth  for  this  event 
produced  a  strong  tradeoff  between  depth  and  yield.  The  bubble  pulse  period  was  consis¬ 
tent  with  the  actual  yield  of  2000  kg.  It  was  found  that  depth  and  yield  could  not  be 
uniquely  constrained  by  the  data  alone.  The  best  depth  and  yield  estimate,  as  might  be 
expected,  as  found  for  the  largest  5000  kg  explosion.  The  results  show  good  agreement 
with  the  know  depths  and  yields,  and  suggest  that  the  cepstrum  inversion  technique  is  a 
useful  method  for  estimating  depth  and  yield  for  underwater  explosions.  For  the  two  lar¬ 
ger  explosions,  we  have  demonstrated  that  the  method  can  estimate  yield  with  an  accu¬ 
racy  of  10  to  15  %. 
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SECTION  4 

UNDERWATER  EXPLOSIONS  IN  THE  BARENTS  SEA  -  THE  OCTOBER  23, 
1997  MURMANSK  AND  AUGUST  12,  2000  KURSK  EVENTS 


This  section  presents  analysis  of  recent  events  in  the  Barents  Sea  presumed  to  be  under¬ 
water  explosions.  The  first  is  a  Russian  underwater  chemical  blast  detonated  near  Mur¬ 
mansk  on  October  23,  1997.  The  second  may  have  been  associated  with  the  Kursk  sub¬ 
marine  sinking  on  August  12,  2000.  The  Murmansk  event  was  described  in  a  note  submit¬ 
ted  to  DTRA  (Baumgardt,  1998)  in  which  it  was  suggested  that  the  event  be  retained  in  a 
database  for  special  event  analysis  of  subsequent  events  in  the  region.  The  August  12, 
2000  Kursk  event,  which  was  also  described  in  another  technical  memo  submitted  to 
DTRA  (Baumgardt  and  Freeman,  2000b),  is  an  example  of  the  kind  of  special  event  that 
would  benefit  from  having  this  kind  of  reference  event  available  for  analysis.  Other  stud¬ 
ies  of  this  event  have  been  published  subsequent  to  ours  (Koper  et  al,  2001;  Gitterman, 
2002),  which  we  will  discuss  later  in  this  section.  In  this  section,  we  first  present  our  re¬ 
analysis  of  the  Murmansk  event.  We  then  discuss  the  August  12,  2000  event  to  show  that 
the  cepstral  characterization  of  the  events  clearly  shows  that  they  are  consistent  with  be¬ 
ing  underwater  explosions  of  different  yield  and  depth. 

4.1  OCTOBER  23,  1997  MURMANSK  EXPLOSION  ANALYSIS. 

The  Murmansk  explosion  of  was  small,  with  a  magnitude  less  than  2.5,  and  was  too  small 
to  be  reported  in  the  REB  of  the  PIDC.  However,  four  seismic  stations,  ARCES,  FINES, 
Spitzbergen,  and  NORES,  recorded  regional  seismic  phases  from  the  event.  Figure  30 
shows  a  record  section  of  the  stations  that  recorded  the  event.  Regional  phases  Pn,  Sn, 
and  Lg  were  recorded  at  each  of  the  stations.  Since  no  location  was  available,  the  event 
was  relocated  using  the  travel  times  for  these  regional  phases.  As  in  the  case  of  the 
Tromso  event  (Section  1.0),  we  used  the  IASPEI91  and  Baltic  Shield  model  to  locate  the 
event.  In  addition,  we  also  obtained  data  from  KEV  to  relocate  the  event.  The  estimated 
locations,  using  different  phases  and  travel-time  curves,  are  shown  in  Figure  3 1 . 
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Figure  30.  Waveforms  recorded  from  the  October  23, 1997  explosion. 


Figure  31.  (Left)  Map  showing  locations  of  the  shot  point  and  the  stations  that  re¬ 
corded  phases  from  the  October  23,  1997  event.  (Right)  Location  esti¬ 
mates  and  error  ellipses  for  the  event. 
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The  locations  in  Figure  31  (right)  all  place  the  event  offshore,  although  some  of  the  error 

{ 

ellipses  overlap  land.  The  best  location,  using  Pn  arrivals  and  the  Baltic  Shield  travel¬ 
time  model,  has  the  smallest  error  ellipse  that  does  not  overlap  land. 

Figure  32  shows  array-stacked  spectra  computed  for  the  regional  phases  shown  in  Figure 
30.  In  each  case,  the  spectra  have  been  corrected  for  the  instrument.  Strong  spectral 
modulations  are  apparent  in  the  spectra  of  the  phases  but  not  in  the  noise,  except  for  those 
spectra  that  are  near  or  below  the  noise  level.  Signal-to-noise  ratios  are  obviously  best  for 
the  ARCES  array,  where  the  entire  modulation  pattern  is  above  the  noise  level. 


Figure  32.  Array  stacked  spectra,  corrected  for  instrument,  for  regional  phases  from 
the  October  23,  1997  Murmansk  event  recorded  at  regional  arrays,  (a) 
ARCES,  (b)  Spitzbergen,  (c)  FINES,  and  (d)  NORES. 
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Because  ARCES  recorded  the  highest  signal  to  noise  ratios,  we  chose  to  model  to  try  to 
model  this  station.  Figure  33  shows  the  spectra  computed  from  the  array-averaged  spectra 
for  ARCES,  using  techniques  discussed  earlier,  and  our  initial  interpretation  of  the  cep- 
stral  features. 


Water  Column 
Reverberations 


Figure  33.  ARCES  cepstra  for  each  phase  computed  from  array-stacked  spectra. 
Cepstra  have  been  shifted  for  display  purposes. 

The  cepstra  are  very  complicated,  with  a  very  evident  negative  trough  at  about  0.32  sec¬ 
onds.  Assuming  this  to  be  the  two-way  acoustic  travel  time  gives  a  depth  of  about  240  m, 
assuming  a  sound  propagation  speed  in  water  of  1480  m/sec,  which  is  close  to  the  known 
bathymetric  water  depths  in  the  region,  which  are  about  200  m.  However,  multiple  posi- 
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tive  peaks  are  observed  which  we  have  interpreted  as  being  caused  by  multiple  bubble 
pulses.  The  question  is  whether  these  assumed  bubble-pulse  delay  times  can  be  accom¬ 
modated  with  the  220  m  depth  in  the  water. 

In  this  study,  we  modeled  the  stack  of  the  phase  cepstra  at  ARCES  using  trial  and  error  to 
check  the  interpretation  of  the  peaks  in  Figure  33.  First,  we  considered  the  possibility' 
that  the  first  peak  at  0.7  seconds  is  a  bubble  pulse  peak.  The  model  construction  for  this  is 
shown  in  Figure  34. 


High  Yield  -  Low  Bottom  Reflection  Coefficient  Model 


Figure  34.  Model  construction  for  a  high-yield  interpretation  of  the  Murmansk  ex¬ 
plosion  cepstra  at  ARCES.  Surface  reflection  coefficient  =  -0.6,  bottom 
reflection  coefficient  =  0.2,  water  depth  =  240  m,  depth  of  explosion  =  240 
m,  yie!d=27,000  kg. 

Figure  35  compares  the  synthetic  cepstrum  for  the  high-yield  model  with  stack  of  all  the 
phase  cepstra  at  ARCES.  The  large  negative  peak  at  about  0.32  seconds  is  matched  well 
by  the  surface  reflection  for  240  m  depth.  However,  to  match  the  large  peak  at  0.7  sec- 
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onds  with  a  bubble  pulse  peak,  the  explosion  yield  must  by  27,000  kg,  which  is  an  order 
of  magnitude  larger  than  the  largest  Dead  Sea  explosion.  If  this  yield  is  right,  the  magni¬ 
tude  of  the  event  should  be  over  4.0.  However,  this  event  had  a  very  low  magnitude  of 
less  than  2.5,  since  it  was  not  reported  by  PIDC  in  the  RJEB.  Thus,  we  are  forced  to  reject 
this  model. 


High  Yield  Cepstral  Model  Fit  to  Data 

Cepstrum 


Quefrency  (sec) 


Figure  35.  Comparison  of  the  high-yield  cepstral  model  with  the  stack  of  all  the 
phase  cepstra  at  ARCES. 

We  note  that  there  is  an  unexplained  peak  at  about  0.2  seconds  which  might  actually  be 
the  bubble  pulse.  Thus,  we  consider  the  possibility  that  this  peak  is  the  bubble  pulse  and 
the  peak  at  0.7  seconds  is  due  to  a  water  reverberation.  Figure  36  shows  a  second,  low 
yield  model.  In  this  model,  we  have  reduced  the  yield  to  200  kg  and  increased  the  bottom 


47 


reflection  coefficient  from  0.2  to  0.6.  Also,  we  consider  two  bubble  pulses. 


Low  Yield  -  High  Bottom  Reflection  Coefficient  Model 


Quefrency  (sec) 


Figure  36.  Model  construction  for  a  high-yield  interpretation  of  the  Murmansk  ex¬ 
plosion  cepstra  at  ARCES.  Surface  reflection  coefficient  =  -0.8,  bottom 
reflection  coefficient  =  0.6,  water  depth  =  240  m,  depth  of  explosion  =  235 
m,  yield=200  kg. 

Figure  37  compares  the  synthetic  cepstrum  with  the  stack  of  all  the  individual  phase  cep¬ 
stra  at  ARCES.  With  this  model,  we  now  explain  the  0.2-second  peak,  as  being  the  pri¬ 
mary  bubble  peak  and  the  peak  at  0.7  seconds  is  now  the  first  bottom  reflection.  How¬ 
ever,  to  get  a  peak  that  large  from  a  bottom  reflection,  we  must  assume  a  bottom  reflec¬ 
tion  coefficient  of  0.6.  This  depth  is  close  to  the  known  water  depth  in  the  region  and  the 
200  kg  yield  estimate  is  more  consistent  with  a  low-magnitude  event. 
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Low  Yield  Cepstral  Model  Fit  to  Data 

Cepstrum 


Figure  37.  Comparison  of  low-yield  synthetic  cepstrum  with  the  stack  of  all  the 
phase  cepstra  at  ARCES. 

One  of  the  problems  we  have  had  in  matching  these  cepstra,  as  noted  earlier,  is  that  a 
small  surface  reflection  coefficient  has  to  be  assumed  to  explain  the  negative  peak.  Nor¬ 
mally,  we  would  expect  this  coefficient  to  be  about  -0.9.  However,  Figure  37  shows  that 
assuming  a  -0.8  surface  reflection  coefficient  produces  a  negative  trough  that  is  much 
larger  than  the  observed  trough.  This  trough  is  produced  by  the  negative  correlation  of  the 
primary  pressure  pulse  in  the  water  and  the  surface  reflection.  The  fact  that  this  peak  is  so 
much  smaller  than  expected  from  the  model  indicates  that  the  two  pulses  are  somewhat 
decorrelated.  This  lack  of  correlation  may  have  been  caused  by  distortion  of  the  reflection 
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off  the  surface,  or  distortion  of  the  signal  as  it  propagates  through  the  water  column.  Our 
model  is  currently  not  able  to  explain  this  decorrelation. 

We  conclude,  based  on  this  analysis,  that  the  Murmansk  event  was  in  fact  a  small  under¬ 
water  explosion,  detonated  on  the  bottom  of  the  Barents  Sea  near  Murmansk,  with  a  yield 
of  about  200  kg. 

4.2  AUGUST  12,  2000  KURSK  EVENT  ANALYSIS. 

In  this  section,  we  interpret  the  seismic  signals  recorded  from  the  Kursk  submarine  acci¬ 
dent  that  occurred  on  August  12,  2000. 

The  Center  for  Monitoring  Research  (CMR)  issued  a  report  (Event  Report,  2000)  describ¬ 
ing  two  events  on  that  date,  a  small  one  (ML=2.2)  and  a  later  larger  one  (mb=3.6).  The 
larger  event  occurred  at  07:30:42  GMT  and  was  described  as  having  characteristics  con¬ 
sistent  with  being  an  underwater  explosion.  In  this  study,  we  further  analyze  the  event 
using  the  cepstral  modeling  and  inversion  algorithm,  discussed  earlier. 

4.2.1  Waveform  Analysis 

Figure  38  shows  on  a  map  the  location  of  the  second  large  Barents  Sea  event,  given  in 
CMR  report  (Event  Report,  2000),  and  the  locations  of  nearby  Scandinavian  arrays  that 
detected  the  event.  Although  the  location  is  based  on  phases  detected  at  ten  stations,  we 
focus  on  the  waveforms  from  the  arrays,  and  primarily,  the  closest  array,  ARCES 
(ARA0).  Figure  39  shows  a  record  section  of  recordings  of  the  event  at  the  regional  ar¬ 
rays  shown  in  Figure  38  and  the  identified  regional  seismic  phases. 


50 


Figure  38.  Map  showing  the  location  of  the  second  and  largest  seismic  event  that  oc¬ 
curred  on  August  12,  2000  at  07:30:42. 


Figure  39.  Record  section  of  single-channel  recordings  from  the  Scandinavian  ar¬ 
rays,  which  recorded  the  August  12,  2000  event  in  the  Barents  Sea.  Each 
waveform  was  bandpass  filtered  in  frequency  bands,  as  indicated. 
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Regional  phases  Pn  and  Sn  could  be  identified  at  all  stations,  with  the  possible  exception 
of  NORSAR  (NAOOO),  although  a  weak  phase  can  be  discerned  at  about  7:38:00  that 
may  be  Sn.  Pg  and  Lg  are  clearly  observed  at  ARCES,  and  Lg  is  observed  at  FINES. 

In  the  remainder  of  this  section,  we  primarily  focus  on  the  ARCES  array  data.  Figure  40 
shows  a  bandpass  filter  analysis  of  the  ARAO  channel  waveform  of  ARCES. 


Figure  40.  Bandpass  filter  analysis  of  the  ARAO  channel  waveform  at  ARCESS. 

The  phase  picks  made  on  the  lower  frequency  bands  are  indicated  with  vertical  lines  that 
also  shows  their  times  at  higher  frequencies.  Strong  Sn  and  Lg  phases  are  evident  at  low 
frequency.  The  Pg  phase  is  very  strong  in  the  2-to-4  and  3-to-5  Hz  band.  However,  at  the 
higher  frequencies,  the  phases  become  much  less  distinct  and  look  more  like  scattered 
waves.  This  result  is  somewhat  different  than  the  October  23,  1997  Murmansk  explosion, 
described  above,  where  we  observed  that  the  explosion  had  distinct  phases  at  the  higher 
frequency.  Another  difference  is  that  the  Murmansk  explosion  had  much  stronger  shear 
waves  relative  to  compressional  waves  at  high  frequency.  Strong  Pn  and  coda  waves  are 
apparent  for  the  August  12,  2000  event  in  Figure  40,  but  the  shear  waves  Sn  and  Lg  are 
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not  as  evident  in  the  coda  at  high  frequency  as  they  were  for  the  Murmansk  event. 

4.2.2  Regional  P/S  Ratios 


The  ratios  of  the  Pn/Sn  and  PgISn  ratios  for  the  event  measured  at  ARCES,  along  with 
measurements  from  earthquakes  and  explosions  in  nearby  regions  also  recorded  at 
ARCES,  are  plotted  as  a  function  of  frequency  in  Figure  41 .  The  earthquakes  occurred  in 
northwestern  Norway  and  in  the  Norwegian  Sea,  the  mine  blasts  on  the  Kola  Peninsula, 
and  the  nuclear  explosions  at  Novaya  Zemlya.  The  marine  explosion  points  (blue  trian¬ 
gles)  are  for  the  presumed  underwater  explosion  or  implosion  associated  with  the  earlier 
Komsomolets  submarine  sinking  in  1989.  A  line  connects  the  points  for  the  August  12, 
2000  event.  All  ratios  have  been  corrected  for  distance. 
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Figure  41.  Pn/Sn  (left)  and  Pg/Sn  ratios  plotted  as  a  function  of  filter  frequency 
band.  The  ratios  have  been  corrected  for  distance. 
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The  connected  points  for  the  August  12,  2000  event  fall  between  the  mine  blast  and  earth¬ 
quake  points  for  the  Pn/Sn  ratios  at  lower  frequencies  but  are  more  explosion-like  at 
higher  frequency.  The  Pg/Sn  ratios  are  clearly  explosion  like  across  the  entire  frequency 
band.  This  trend  can  also  be  clearly  seen  in  Fig,  where  the  Pg  and  Sn  phases  are  both 
quite  strong  at  low  frequency,  and  Pg  is  very  impulsive  and  larger  than  Sn  in  the  0.5-2. 5 
and  2-4  Hz  bands. 

Overall,  the  Pn/Sn  and  Pg/Sn  ratios  are  consistent  with  the  event  being  an  underwater 
explosion. 

4.2.3  Spectral/Cepstral  Analysis 

Figure  42  shows  the  array-averaged  spectra  for  each  regional  phase  and  noise. 


Figure  42.  Array-averaged  spectra  for  each  regional  phase  from  the  12  August  2000 
event  and  the  background  noise  to  Pn  recorded  on  the  vertical  component 
channels  of  ARCES. 
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These  spectra  are  computed  by  first  making  spectral  estimates  of  each  of  the  four  phases, 
Pn,  Pg,  Sn,  and  Lg,  on  each  vertical  component  channel  of  the  array  and  then  averaging 
the  spectra  across  all  array  elements.  The  noise  spectra  ahead  of  the  Pn  phase  on  each 
channel  are  computed  in  the  same  way.  The  spectra  have  not  been  corrected  for  instru¬ 
ment  response.  Thus,  the  spectra  go  to  low  values  at  low  frequency,  because  of  the  low- 
frequency  instrument  response  roll  off,  and  at  about  18  Hz  on  the  high  end  due  to  the 
anti-aliasing  filter.  The  most  notable  feature  of  the  spectra  is  the  strong  modulations  or 
“scalloping”  of  the  spectra  that  appear  in  all  phases.  The  black  vertical  dashed  lines  indi¬ 
cate  the  nulls  of  essentially  two  interfering  scalloping  patterns.  The  broader  scalloping, 
indicated  by  the  black  vertical  lines,  have  nulls  separated  by  about  5  Hz.  This  broad  scal¬ 
loping  results  for  spectral  interference  of  reflections  in  the  water  column.  The  first  peak 
of  this  broad  scalloping  has  another  scalloping  pattern  superimposed,  as  indicated  by  the 
red  dashed  lines,  due  to  the  bubble  pulse. 

Figure  43  shows  the  cosine  or  signed  cepstra  computed  from  the  spectra  in  Figure  42. 


Figure  43.  Cepstra  computed  from  the  spectra  in  Figure  42.  The  cepstra  for  the 
phases  and  noise  have  been  shifted  for  display  purposes. 
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The  cosine  cepstra  are  computed  by  methods  discussed  above.  The  red  dashed  line  indi¬ 
cates  the  primary  positive  cepstral  peak  in  the  phase  cepstra.  Note  that  this  peak  does  not 
appear  in  the  noise.  We  interpret  this  peak  as  resulting  from  the  bubble-pulse  scalloping, 
as  indicated  in  Figure  42,  and  has  a  delay  time  of  about  0.75  seconds.  The  black  dashed 
line  indicates  a  negative  peak  at  about  0.2  seconds  that  also  does  not  appear  the  noise.  We 
interpret  this  peak  as  resulting  from  the  broad  spectral  scalloping  in  Figure  42  caused  by 
the  surface  reflection,  which  would  have  a  negative  polarity  relative  to  the  primary  phase. 
Thus,  we  conclude  that  these  spectra  and  cepstra  are  consistent  with  the  characteristics  of 
underwater  explosions. 

4.2.4  Cepstral  Model  and  Inversion 

To  further  confirm  that  this  event  was  an  underwater  explosion,  and  to  estimate  its  depth 
and  yield,  we  attempt  to  fit  model  cepstra  for  explosions  to  the  observed  cepstra. 

For  modeling  purposes,  we  again  choose  to  model  the  “stacked”  cepstrum,  that  is,  the  av¬ 
erage  of  each  of  the  individual  phase  cepstra.  The  stacked  cepstrum  provides  some  im¬ 
provement  in  signal  strength  and  reduces  features  that  do  not  appear  in  all  the  cepstra. 
Because  the  model  cepstra  are  phase- independent,  the  stacked  cepstrum  serves  as  the  best 
overall  representation  of  the  observed  cepstra. 

In  Figure  44,  both  the  spectra  and  cepstra  are  replotted  with  the  stack  spectrum  and  cep¬ 
strum  plotted  as  the  green  trace.  In  this  case,  the  spectra  and  cepstra  of  the  individual 
phases  are  very  similar  so  that  the  stacked  spectra  and  cepstra  are  very  similar  to  the  indi¬ 
vidual  ones.  From  now  on,  we  use  the  stacked  cepstrum  as  the  observed  cepstrum  that  we 
attempt  to  fit  with  a  model  cepstrum. 
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Figure  44.  (Left)  Noise,  phase  and  stacked  spectra;  (right)  Phase,  noise  and  stacked 
cepstra.  The  green  traces  are  the  stacked  spectra  and  cepstra. 

Also,  we  are  only  interested  in  modeling  features  that  appear  in  all  of  the  individual  phase 
cepstra,  and  not  in  the  noise  cepstra.  Generally,  we  regard  any  feature  that  appears  in  both 
the  signal  and  noise  cepstra  to  be  due  to  effects  of  noise  or  processing  artifact.  For  exam¬ 
ple,  in  Figure  43,  we  point  out  apeak  that  appears  in  the  cepstra  just  ahead  of  the  bubble- 
pulse  peak.  This  peak  may  also  be  in  the  noise  cepstra,  so  we  choose  to  ignore  it  when 
fitting  model  and  signal  cepstra.  In  order  to  try  to  reduce  noise  effects,  we  have  experi¬ 
mented  with  two  noise  subtraction  methods.  In  the  first  method,  we  subtract  noise  spectra 
from  each  of  the  individual  phase  spectra  before  the  cepstrum  computation.  In  the  second 
method,  we  subtract  the  noise  cepstrum  from  the  signal  cepstrum.  In  this  study,  we  show 
the  results  of  fitting  modeled  cepstra  to  observed  cepstra  corrected  by  spectral  noise  sub¬ 
traction.  The  same  results  were  obtained  with  cepstral  noise  subtraction. 

We  use  the  exhaustive  grid-search  method  for  matching  synthetic  and  observed  cepstra 
and  use  the  one  minus  the  LI  norm  ( 1  -LI )  as  the  match  statistic,  where  the  LI  nonn  is  the 
averaged  absolute  value  of  the  difference  between  the  observed  and  synthetic  cepstra. 
plot.  Figure  45  shows  the  result  of  this  analysis  on  the  spectral-noise  corrected  cepstra. 
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Figure  45.  Contour  plot  of  1-L1  as  a  function  of  depth  and  yield  for  matches  of  syn¬ 
thetic  cepstra  to  observed  cepstra  corrected  by  using  noise  spectra. 

This  plot  shows  that  the  solution  is  non-unique  in  that  there  are  two  peaks,  one  at  low 
yield,  about  4500  kg,  and  the  second  at  a  considerably  higher  yield  of  over  10000  kg.  The 
higher-yield  peak  has  a  slightly  higher  1-L1  statistic  which  is  why  it  was  picked  in  Figure 
45.  Also,  the  depth  of  the  lower- yield  peak  is  about  90  m  compared  with  a  depth  of  about 
130  m  for  the  higher- yield  peak.  Both  of  these  depths  are  reasonably  consistent  with  the 
known  bathymetry  in  the  region,  where  the  water  depths  are  on  the  order  of  250  m. 


Figure  46  shows  cross  sections  at  constant  depth  through  these  two  peaks.  The  width  of 
these  peaks  gives  some  measure  of  the  yield-estimation  uncertainty  of  this  method. 
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Figure  46.  Constant-depth  cross  sections  through  the  two  peaks  shown  in  Figure  45. 

Both  peaks  are  significant  and  their  1  -LI  values  are  nearly  the  same,  or  about  0.97. 
Therefore,  these  plots  alone  are  not  sufficient  alone  to  determine  which  model  is  best. 
Figure  47  shows  a  comparison  of  the  synthetic  and  observed  stacked  cepstra  for  these  two 
peaks. 


Figure  47.  (Left)  Comparison  of  synthetic  and  observed  cepstra  for  the  low-yield 

peak  in  Figure  45.  (right)  Comparison  of  synthetic  and  observed  cepstra 
for  the  high-yield  peak  in  Figure  45. 

Both  models  match  the  bubble-pulse  cepstral  peak  amplitude  and  delay  time  reasonably 
well.  Looking  at  the  spectra  in  Figure  42,  there  to  be  two  broad  modulations  due  to  rever¬ 
berations  with  peak  separations  of  4  Hz  and  6  Hz.  The  former  produces  a  cepstral  trough 
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at  0.25  Hz  and  the  latter  produces  one  at  0.16  seconds.  Both  of  the  troughs  appear  in  the 
cepstrums  in  Figure  43.  The  depth  of  90  m  appears  to  be  matching  the  0.16-second  trough 
on  the  stack  cepstrum,  but  not  the  0.25  trough,  whereas  the  130  m  depth  is  matching  the 
0.25  second  trough  by  not  the  0.16.  The  known  water  depths  in  the  region  are  about  100 
m.  Thus,  we  would  favor  the  90  m  depth  with  the  lower  yield. 

One  possible  explanation  for  the  two  sets  of  modulations  is  that  the  4  Hz  modulation  ex¬ 
tends  from  about  12  Hz  to  16  Hz,  and  the  antialias  filter  rolloff  begins  at  about  16  Hz.  It 
could  that  the  apparent  4  Hz  modulation  is  due  to  the  filter  cutting  off  the  high-frequency 
lobe  of  the  modulation.  Thus,  this  modulation  frequency  may  be  lower  than  it  should  be, 
producing  the  0.25  Hz  peak  in  the  cepstrum. 

Another  possibility  is  that  the  surface  modulation  may  be  distorted  by  the  proximity  of 
the  explosion  to  other  reflection  surfaces,  such  as  the  submarine  itself,  which  could  dis¬ 
tort  the  spectral  modulation.  Other  spectral  peaks  are  evident  in  the  spectra  in  Figure  42  at 
frequencies  above  16  Hz,  which  may  be  causing  some  of  this  distortion.  So,  we  believe 
that  matching  the  0.25-second  cepstral  peak  may  not  be  that  important  since  it  may  be 
distorted.  The  90  m  depth  model  seems  to  match  the  very  strong,  low-quefrency  peaks 
more  accurately.  Thus,  we  believe  that  this  match  is  more  likely  the  correct  model. 

4.3  CONCLUSIONS. 

We  conclude  overall,  based  on  our  analysis,  that  both  events  in  the  Barents  Sea  that  oc¬ 
curred  on  October  23,  1997  and  August  12,  2000  were  underwater  blasts.  For  the  Mur¬ 
mansk  event,  we  used  a  trial  and  error  matching  method  to  obtain  synthetic  cepstrum  that 
explains  the  essential  characteristics  of  the  observed  cepstrum.  The  lower  yield  of  the 
event,  which  was  about  200  kg,  with  depth  of  235  to  240  m  is  consistent  with  the  small 
magnitude  of  the  event  and  the  known  bathymetry  in  the  region. 

The  second  event  on  August  1 2,  2000  has  signal  characteristics  consistent  with  the  event 
being  an  underwater  explosion,  similar  to  those  observed  for  the  Murmansk  event.  Using 
the  cepstral  modeling  and  matching  method,  we  obtain  a  best-matching  model  that  has  a 
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yield  of  10000  kg,  or  about  10  tons,  and  a  depth  of  130  m.  However,  we  also  obtained  a 
low  yield  fit  of  4300  kg,  or  4.73  tons,  at  a  depth  of  90  m.  This  lower-yield  result  seems  to 
be  more  consistent  with  the  expected  yield  of  a  torpedo  as  well  as  the  known  depth  of  wa¬ 
ter  in  the  region,  which  is  about  10  m. 

News  reports  of  sonar  data  indicated  that  the  submarine  may  have  been  sunk  by  a  very 
large  explosion,  perhaps  due  to  an  accidental  torpedo  explosion.  It  is  possible  that  the 
explosion  itself  happened  inside  the  submarine  before  it  sank  to  the  bottom,  or  the  subma¬ 
rine  may  have  rested  somewhat  vertically  on  the  ocean  bottom.  Thus,  the  explosion  may 
have  occurred  on  or  inside  the  submarine  but  elevated  above  the  water  bottom.  Our  re¬ 
sults  indicate  that  if  the  90  m  depth  is  right,  the  explosion  occurred  about  10  m  above  the 
surface. 

Our  estimates  of  yield  essentially  agree  with  the  results  of  Koper  et  al  (2001).  However, 
Gitterman  (2002),  using  some  of  the  same  data  we  have  analyzed,  recently  published  a 
lower  yield  of  1500  and  1650  kg  for  depths  of  61  and  64  m.  He  states  that  the  bubble  fre¬ 
quency  was  1 .45  Hz,  using  ARCES  data,  which  is  close  to  our  estimates.  However,  he 
uses  the  depths  of  61  and  64  m  to  estimate  the  yields.  After  repeated  readings  of  this  pa¬ 
per,  we  cannot  determine  how  he  arrived  at  these  depths.  Apparently,  they  were  not  de¬ 
termined  from  data  since  he  states  that  he  was  not  able  to  observe  any  evidence  of  spec¬ 
tral  modulations  due  to  the  water  reverberations.  However,  we  have  shown  clear  evidence 
in  this  report  of  spectral  modulations  due  to  reverberations.  Thus,  we  argue  that  Gitter¬ 
man  (2002)  has  underestimated  the  yields  because  his  “assumed”  detonation  depths  of  6 1 
to  64  m  are  too  low  and  are  not  supported  by  the  data. 
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SECTION  5 

OVERALL  CONCLUSIONS  AND  RECOMMENDATIONS 


This  study  has  shown  that  seismic  recordings  of  underwater  blasts  can  be  useful  for  the 
estimation  of  yield  and  water  depths  of  the  explosions.  We  have  shown  that  the  model 
developed  by  Baumgardt  and  Der  (1998)  effectively  matches  the  predominant  features  of 
cepstra  of  underwater  explosions,  and  simple  fitting  methods,  using  search  techniques, 
can  “invert”  observed  cepstra  to  infer  the  depth  and  yield  of  the  explosion.  The  cepstral 
analysis  technique  provides  a  method  for  unambiguously  separating  out  effects  of  bubble 
pulses  and  reverberation  spectral  modulations,  and  is  a  simpler  technique  that  trying  ana¬ 
lyze  harmonic  peaks  (e.g.,  Gitterman,  2002).  Thus,  the  method  can  easily  be  automated 
for  routine  processing  of  seismic  data. 

The  method  described  in  this  report  was  prototyped  in  Matlab©  and  is  described  more  in 
a  users  manual  that  can  be  found  in  the  Appendix.  We  recommend  that  this  technique  be 
incorporated  as  part  of  screening  of  events  that  are  located  off  shore,  including  those  that 
occur  close  to  coasts  where  the  error  ellipses  overlap  land  areas.  If  an  event  that  occurs 
offshore  can  be  fit  with  an  underwater  explosion  cepstral  model,  then  it  could  be  so  iden¬ 
tified.  Our  analyses  of  many  offshore  events  near  Scandinavia  in  our  first  report 
(Baumgardt,  1 999)  showed  that  many  offshore  events  are  underwater  blasts  and  the 
automated  screening  methods  being  used  by  the  IDC  and  other  national  data  centers 
should  properly  identify  such  events. 

More  research  is  required  to  refine  the  matching  techniques  and  to  estimate  confidence  of 
match.  Studies  of  how  this  method  works  on  earthquake  spectra  should  also  be  made.  In 
our  Annual  Report  (Baumgardt,  1999),  we  showed  some  data  that  indicated  that  undersea 
earthquakes  could  produce  spectral  modulations  due  to  water  reverberations.  Discrimina¬ 
tion  studies  might  identify  these  events,  however  our  analysis  of  discriminants  for  under¬ 
water  blasts  suggest  that  many  of  them  are  very  earthquake  like.  Thus,  the  cepstral  match¬ 
ing  method  may  be  needed  for  discrimination.  Thus,  methods  for  how  to  evaluate 
matches  to  earthquakes,  or  lack  of  good  matches,  must  be  developed. 
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Of  course,  our  studies  have  been  limited  to  chemical  blasts.  Generally,  the  yields  of  these 
events  are  very  low,  compared  to  those  for  even  small  nuclear  blasts.  The  task  of  monitor¬ 
ing  for  underwater  blasts  currently  resides  with  the  hydroacoustic  networks,  and  the 
methods  described  in  this  report,  with  some  modification,  could  also  be  applied  to  hy¬ 
droacoustic  data.  However,  we  recommend  that  seismic  monitoring  systems  also  have  the 
capability  of  characterizing  underwater  blasts,  since  explosions  may  occur  in  inlets  or 
inland  seas  which  may  cause  hydroacoustic  signals  to  be  blocked.  Also,  as  we  argued  in 
the  introduction,  seismic  signals  propagate  more  directly  from  the  source  than  do  hy¬ 
droacoustic  signals.  Because  underwater  blasts  will  likely  produce  sizable  seismic  sig¬ 
nals,  techniques  such  as  this  should  be  made  available  to  the  seismic  monitoring  systems 
as  well  as  hydroacoustic  systems. 
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A 1  GRAPHICAL  USER  INTERFACE  (GUI) 

The  cepstral  analysis  modeling  and  inversion  tool  for  underwater  explosion  allows  the 
user  to  create  and  save  cepstra  files  from  selected  data,  and  to  run  an  inversion  program 
on  previously  processed  cepstra  files.  The  inversion  determines  a  best  match,  using 
various  selectable  comparison  and  optimization  methods  between  the  data  cepstrum  and 
various  theoretical  cepstra.  For  both  of  these  options,  input  parameters  must  be  specified. 

Al.l  MODELING  AND  INVERSION  PROGRAM. 

The  user  interface  for  the  Cepstral  Modeling  and  Inversion  Tool  is  Windows-based,  and 
was  designed  using  Matlab. 

To  start  the  program,  double-click  the  Cepstral  Analysis  icon.  The  Cepstral  Modeling 
and  Inversion  Program  Window  appears  on  the  screen  as  shown  in  Figure  A-l.  The 
Cepstral  Inversion  Program  may  also  be  started  by  opening  matlab,  and  typing  cep_ana 
on  the  command  line. 

Note  that  the  location  of  the  cepstral  inversion  program  needs  to  be  added  to  the  Matlab 
path  for  the  program  to  run  from  any  directory  in  Matlab.  The  Cepstrum  folder, 
containing  the  cepstral  inversion  algorithm,  should  be  placed  in  the 
/matlab/ toolbox/Cepstrum  directory.  If  the  . .  ./Cepstrum  folder  is  not  added  to  the  matlab 
path,  the  user  must  be  in  the  directory  of  the  cepstral  inversion  program  for  the  command 
cep_ana  to  run  the  inversion  program. 
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Figure  A-l.Cepstral  Analysis  Modeling  and  Inversion  Tool  for  Underwater  Explo¬ 
sions  Window 

The  Cepstral  Analysis  Modeling  and  Inversion  Program  window  has  two  menus: 
Parameters  and  Select. 

Al.1.1  Parameters  Menu 

The  Parameters  menu  has  one  option:  Stations  and  Channels.  A  detailed  description 
of  this  option  is  given  with  Figure  A-2.  The  Stations  and  Channels  option  is  accessed 
by  a  pull-down  option  from  the  Parameters  menu.  Selections  in  the  Stations  and 
Channels  window  must  be  made  in  order  to  compute  and  save  cepstra  files  from  data. 


Parameters 


Stations  and  Channels 


Figure  A-2.  Parameters  Menu  and  Options 

Stations  and  Channels:  Displays  a  window  containing  inputs  for  the  waveform  (*.w) 
file  and  inprefix  of  the  additional  files  needed  for  processing  the  data.  The  files  needed 
are: 

*.w  Waveform  File 
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*. affiliation  Network  station  affiliations 

*. arrival  Summary  information  on  a  seismic  arrival 

*. assoc  Data  associating  arrivals  with  origins 

*. origin  ata  on  event  location  and  confidence  bounds 

*.wfdisc  aveform  file  header  and  descriptive  information 

*.wftag  Waveform  mapping  file 

Once  these  two  inputs  are  chosen,  the  stations  and  channels  for  the  particular  waveform 
file  chosen  are  displayed  (section  A1.2). 

Al.1.2  Select  Menu 

The  Select  menu  (Figure  A-3)  has  two  options:  Data  Cepstrum  and  Inversion.  A 
detailed  description  of  these  two  options  is  given  in  sections  A1 .3  and  A1 .4. 


D  ala  Cepstrum 
Inversion 


Figure  A-3.  Select  Menu  and  Options 


A1.2  STATIONS  AND  CHANNELS. 

The  Stations  and  Channels  GUI  window  is  shown  in  Figure  A-5. 

Waveform  Files...:  Clicking  on  this  input  button  opens  a  browse  window  which  filters 
for  *.w  files  (Figure  A-4).  It  is  possible  to  browse,  and  also  to  select  a  different  folder 
containing  waveform  (.w)  files  by  changing  the  Look  in  folder  in  the  browse  window. 
Select  the  desired  *.w  file  to  process,  either  by  selecting  the  file  in  the  browse  window, 
and  selecting  Open,  or  by  double  clicking  on  the  desired  *.w  file  in  the  browse  window. 
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Figure  A-4.  Choose  .w  File  to  Analyze  Browse  Window 

Prefix  for  Files:  This  is  the  prefix  of  the  files  needed,  in  addition  to  the  *.w  file  for 
processing  the  data  (section  A1 .1 . 1).  Often,  the  prefixes  for  the  *.w  file  and  the  other 
files  are  identical.  Therefore,  this  edit  box  will  fill  in  automatically  with  the  prefix  of  the 
*.w  file  selected.  However,  when  the  prefixes  for  the  additional  files  and  the  .w  file 
differ,  the  Prefix  for  Files  input  box  is  edited  to  reflect  the  inprefix  name  of  the 
additional  files. 

Stations: 

Select  All:  Checks  all  of  the  checkboxes  for  stations  listed  in  the  window. 

Clear  All:  Clears  all  of  the  checkboxes  for  the  stations  list  in  the  window. 

The  station  names  for  the  chosen  .w  file(s)  are  listed  (EIL,  MLR,  MRNI,  Figure  A-5)  in  a 
row.  The  elements  of  the  respective  stations  are  listed  in  up  to  two  columns  under  the 
station  name.  The  channels  for  each  of  the  stations  are  listed  under  the  station  elements 
(be,  bn,  bz.  Figure  A-5).  At  least  one  element  and  one  channel  for  one  of  the  stations 
must  be  chosen  in  order  to  compute  the  data  cepstrum. 

Stations  and  Channels:  Press  this  pushbutton  after  selecting  waveform  and  inprefix  for 
files.  Stations  and  Channels  will  be  displayed. 
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Save  Selection:  Saves  the  stations  and  channels  to  be  processed.  This  also  saves  the 
selected  *.w  file  and  inprefix  selected  as  defaults  in  the  file,  sta_cha.mat  in  the 
/matlab/toolbox/Cepstrum/MATFELES  directory.  Subsequently,  they  will  be  used  as 
defaults. 

Close:  Closes  the  Stations  and  Channels  window. 
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Figure  A-5.  Stations  and  Channels  Window 


A  1.3  PROCESS  CEPSTRUM. 

The  GUI  window  for  Process  Cepstrum  is  shown  in  Figure  A-5.  The  Process  Cepstrum 
window  is  accessed  by  selecting  Data  Cepstrum  from  the  pull-down  Select  menu  in  the 
Modeling  and  Inversion  window  (section  Al.1.2). 
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Utilities:  The  pull  down  menu  contains  the  option  to  Save  Defaults.  Selecting  Save 
Defaults  saves  the  current  inputs  in  the  Process  Cepstrum  window  for  future  use  as 
defaults  in  a  default_data.mat  file  in  the  /matlab/toolbox/Cepstrum/MATFILES  directory. 

File  Processing:  The  corresponding  box  displays  the  *.w  file  currently  processing. 
Clicking  on  the  File  Processing...  pushbutton  opens  a  browse  window.  Existing  *.w 
files  may  be  browsed;  however,  if  a  different  *.w  file  is  to  be  processed  from  the  one 
displayed,  the  change  of  .w  file  must  be  made  in  the  Stations  and  Channels  window 
(section  A1 .2).  This  window  will  also  display  the  corresponding  stations  and  channels  of 
the  .w  file  for  the  user  to  select. 

Inprefix:  The  inprefix  for  the  additional  files,  besides  the  *.w  file  used  to  process  data  is 
displayed.  The  edit  box  contains  the  prefix  for  files  specified  in  the  Inprefix  edit  box  in 
the  Stations  and  Channels  GUI. 

Window:  Specifies  type  of  spectral  smoothing.  The  windowing  options  to  apply  to  the 
data  are  Hanning,  Kaiser,  or  Parzen. 

Length  (s):  The  edit  box  specifies  the  length  of  the  spectral  window,  in  seconds. 

Freq  From:  The  edit  box  specifies  the  low  frequency,  in  hertz  (Hz),  to  filter  the  signal. 
This  is  the  lower  end  of  the  bandwidth  to  compute  the  cepstrum 

To:  The  edit  box  specifies  the  high  frequency  to  filter  to,  in  Hz.  This  is  the  higher  end  of 
the  bandwidth  to  compute  the  cepstrum. 

Subtract  Noise:  If  checked,  noise  is  subtracted  from  each  phase.  Noise  is  defined  as  the 
time  before  the  first  phase  available  for  the  station  and  channel.  The  noise  is  processed  in 
the  same  way  as  the  phase  data. 

Remove  Instr.  Resp.:  If  checked,  the  instrument  response  is  removed  by  dividing  out  the 
response.  The  instrument  response  is  read  in  from  a  file. 

Plot  Raw  Spectra:  If  checked,  the  spectra  of  each  phase  and  channel,  without 
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polynomial  subtraction  or  application  of  frequency  limits  specified  in  the  Freq  From: 
and  To:  edit  boxes,  are  displayed,  in  a  separate  figure  for  each  station. 

Plot  Spectra:  If  checked,  the  processed  spectra  for  each  phase  and  channel  are  displayed 
in  a  separate  figure  for  each  station. 

Polynomial  Fit  Degree:  Removes  the  trend  from  the  spectra.  The  contents  of  the  edit 
box  specify  which  degree  polynomial  to  fit  to  the  spectra.  This  trend  is  then  subtracted 
from  the  data. 

Cepstrum: 

Stacked  Spectrum:  Compute  the  data  cepstrum  by  stacking  the  processed  spectra,  then 
taking  the  cepstrum  of  the  stacked  spectra.  The  cepstra  files  computed  are  saved  in  the 
folder  specified  in  the  Save  Files  input  box,  with  the  extension  .cp. 

Stacked  Cepstrum:  Compute  the  data  cepstrum  by  computing  the  cepstra  of  the 
individual  spectra,  then  stacking  the  cepstra.  The  cepstra  files  computed  are  saved  in  the 
folder  specified  in  the  Save  Files  input  box,  with  the  extension  .sp. 

Cepstrum  Len.  (s):  Length  of  time  to  compute  data  cepstrum  (in  seconds). 

Save  Files:  Location  to  save  cepstra  files,  and  also  variable  file  containing  inputs  to 
creating  the  data  cepstrum  needed  for  inversion  of  data  cepstrum  and  model  cepstrum. 

The  default  for  this  file  location  is  the  . .  ./cepstra_files  folder  in  the  folder  which  contains 
the  data  for  processing.  The  file  containing  the  parameters  is  saved  as  an  ‘inprefix’. mat 
file  in  the  same  folder. 

Process  Cepstrum:  When  selected,  the  variables  specified  in  the  Process  Cepstrum 
window  are  used  to  compute  the  data  cepstrum.  The  calculated  cepstra  are  saved,  both  as 
averaged  individual  phases,  and  as  a  stack.  The  individual  phases  are  saved  as 
inprefix. phase_name.cp  (or  inprefix. phase_name.sp  is  Stacked  Spectrum  was  selected  in 
Process  Cepstrum  GUI).  The  stack  is  saved  as  inprefix. ST. cp  (or  inprefix. ST. sp  if 
Stacked  Spectrum  was  selected  in  the  Process  Cepstrum  GUI).  These  files  are  saved 
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in  the  folder  specified  in  the  Save  Files  input  box.  Also  saved  is  the  file  containing  the 
input  parameters  needed  for  the  inversion,  as  ‘inprefix’.mat. 

Close:  When  selected,  the  Process  Cepstrum  window  is  closed. 


Figure  A-6.  Process  Cepstrum  Window 


A  1.4  INVERSION. 

The  GUI  window  for  Inversion  is  shown  in  Figure  A-7.  The  Inversion  window  shows  a 
saved  file  from  a  previously  saved  processed  cepstrum.  Note  that  the  name  (and  data)  is 
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different  than  that  previously  shown  for  explaining  processing  the  cepstrum  (section 
A1.3).  The  Inversion  window  is  accessed  by  selecting  Inversion  from  the  pull  down 
menu  in  the  Modeling  and  Inversion  window  (section  Al.1.2). 

Utilities:  The  pull  down  menu  contains  the  option  to  Save  Defaults.  This  saves  the 
current  inputs  in  the  Inversion  window  for  future  use  as  an  inv.mat  file  in  the 
/matlab/toolbox/Cepstrum/MATFILES  directory. 

Saved  Files...:  When  selected,  displays  a  browse  window  on  screen  similar  to  that 
shown  in  Figure  A-4.  Clicking  on  this  input  button  opens  a  browse  window  which  filters 
for  *ST.cp  files,  similar  to  Figure  A-4.  It  is  possible  to  browse  and  select  a  different 
folder  containing  cepstra  (*ST.cp)  files  by  changing  the  Look  in  folder  in  the  browse 
window.  It  is  also  possible  to  select  cepstra  files,  which  were  calculated  using  the 
Stacked  Spectrum  option  from  Process  Cepstrum  GUI  window.  These  files  have  a 
Asp  suffix.  Select  the  desired  Acp  (or  Asp)  file  to  process,  either  by  selecting  the  file  in 
the  browse  window,  and  selecting  Open,  or  by  double  clicking  on  the  desired  *.cp  (or 
Asp)  file  in  the  browse  window. 

Inprefix:  Specifies  the  inprefix  to  the  saved  parameters  file  (variables  such  as  sampling 
rate)  saved  to  the  folder  specified  in  the  Save  Files  edit  box  when  the  cepstra  were 
processed  using  the  Process  Cepstrum  GUI.  The  edit  box  defaults  to  the  inprefix  of  the 
saved  cepstrum  file  to  process. 

Comparison: 

Cross  Corr:  If  selected,  the  cross  correlation  between  the  data  cepstra  and  theoretical 
cepstrum  is  computed. 

. llnorm:  If  selected,  llnorm  between  the  data  cepstrum  and  theoretical  cepstra  is 

computed 

. 12norm:  If  selected,  12norm  between  the  data  cepstrum  and  theoretical  cepstra  is 

computed 
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Optimization: 


Range  of  Values:  Exhaustive  Search  Technique  that  steps  through  user  defined  yield  and 
depth  space.  At  each  grid  point  of  yield  and  depth,  the  theoretical  cepstrum  is  computed, 
and  compared  to  the  data  cepstrum  using  the  Comparison  method  chosen. 

Downhill  Simplex:  Optimal  Search  Technique  (originally  developed  by  Nelder  and 
Mead  (1965))  that  does  a  multidimensional  search  of  a  parameter  space  to  minimize  a 
function.  In  this  case,  to  minimizes  the  difference  between  the  data  cepstrum  and 
theoretical  cepstra  using  the  Comparison  method  chosen.  This  method  finds  a  local 
minimum 

Sim.  Annealing:  Optimal  Search  Technique  that  finds  a  global  minimum  between  the 
data  cepstrum  and  theoretical  cepstra,  using  the  comparison  choice  selected  under 

Comparison. 

Temps:  The  edit  box  contains  the  number  of  temperatures  to  be  used  for  the  Simulated 
Annealing  search  technique.  This  box  is  editable  only  if  Sim.  Annealing  is  chosen  as 
Optimization  technique. 

Trials:  The  edit  box  contains  the  number  of  trials  to  be  used  for  the  Simulated 
Annealing  search  technique.  This  box  is  editable  only  if  Sim.  Annealing  is  chosen  as 
Optimization  technique. 

Yield(kg):  Yield  in  Kilograms(kg)  for  theoretical  model.  For  Range  of  Values 
Optimization  method,  this  is  the  beginning  range  of  yield.  For  Downhill  Simplex  and 
Sim.  Annealing  Optimization  methods,  Yield(kg)  is  the  starting  point  for  the 
optimization  method. 

To:  For  Range  of  Values  Optimization  method,  this  is  the  high  end  of  the  yield  range, 
in  kg,  for  the  theoretical  cepstrum.  This  option  is  not  visible  if  Range  of  Values  is  not 
the  chosen  Optimization  method. 
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Step:  The  edit  box  signifies  the  successive  steps  from  Yield(kg)  and  the  value  in  edit  box 
To:  (high  yield  value)  to  take.  This  option  is  only  available  for  the  Range  of  Values 
Optimization  method. 

Depth(m):  Depth  in  meters(m)  for  theoretical  model.  For  Range  of  Values 
Optimization  method,  this  is  the  beginning  range  of  depth.  For  Downhill  Simplex  and 
Sim.  Annealing  Optimization  methods,  Depth(m)  is  the  starting  point,  for  the 
optimization  method. 

To:  For  Range  of  Values  Optimization  method,  this  is  the  high  end  of  the  depth  range, 
in  m,  for  the  theoretical  cepstrum.  This  option  is  not  visible  if  Range  of  Values  is  not  the 
chosen  Optimization  method. 

Step:  The  edit  box  signifies  the  successive  steps  from  Depth(kg)  and  the  value  in  edit 
box  To:  (high  depth  value)  to  take.  This  option  is  only  available  for  the  Range  of  Values 
Optimization  method. 

Water  Depth(m):  The  depth  of  the  water,  in  m,  at  the  location  of  the  presumed 
explosion. 

Sound  Speed  (m/s):  Speed  of  sound  in  water,  in  units  of  (m/s). 

Depth  &  Yield  from  Data  Cepstrum:  If  this  is  checked,  then  an  estimate  of  the 
beginning  depth  and  yield  values  for  the  theoretical  cepstrum  are  calculated  from  user 
input  on  the  graph  of  the  data  cepstrum.  When  Process  Inversion  is  pressed,  a  figure 
containing  the  stacks  of  the  individual  phases  of  the  data  cepstrum,  as  well  as  a  stacked 
cepstrum  over  all  phases  is  shown.  With  the  mouse,  hold  cursor  over  the  surface 
reflection  and  click  once.  Next,  hold  cursor  over  the  first  bubble  pulse  and  click  once. 
The  algorithm  computes  a  starting  depth  and  yield  from  these  inputs. 

Plot  Theoretical  Model:  If  this  is  checked,  then  the  theoretical  model  of  the  bubble 
pulse  wavelet,  water  column  reflection  wavelet,  bubble  pulse  and  reflection  wavelet, 
minimum  phase  wavelet,  and  cosine  cepstrum  will  be  displayed,  for  the  theoretical  model 
which  is  the  best  match  to  the  data  cepstrum  by  the  chosen  Comparison  model. 
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Comparison  to  Theoretical(s):  The  edit  box  specifies  the  time,  in  seconds  to  compare 
the  data  cepstrum  to  the  theoretical  cepstrum. 

Select  Theoretical  Parameters:  If  this  is  selected,  a  Theoretical  Parameters  window  is 
shown  (see  section  0).  This  window  allows  the  user  to  change  theoretical  cepstrum 
parameters. 

Process  Inversion:  When  selected,  the  algorithm  utilizes  the  input  parameters  specified 
in  the  Inversion  window  to  process  the  inversion  of  the  data  cepstrum  with  theoretical 
cepstra. 

Close:  When  selected,  the  Inversion  window  is  closed. 
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Figure  A-7.  Inversion  Window 


A1.5  Theoretical  Cepstrurn  Parameters. 

The  GUI  window  for  Theoretical  Cepstrurn  Parameters  is  shown  in  Figure  A-8.  The 
Theoretical  Cepstrurn  Parameters  window  is  accessed  by  pressing  Select  Theoretical 
Parameters  in  the  Inversion  window.  Default  inputs  are  initially  displayed,  and  used  as 
parameters. 
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Surface  Reflection:  Reflection  from  the  surface  of  the  water  (may  range  in  value  from  - 
1.0  to  0.0). 


Bottom  Reflection:  Reflection  from  the  bottom  of  the  sea  floor  (may  range  in  value  from 
0.0  to  1.0) 

Decay  Constant  for  Bubble  Pulse:  Multiplicative  amplitude  factor,  which  was  set  to 
0.4,  since  this  value  was  required  to  produce  the  spectral  modulation  observed 
(Baumgardt  and  Der,  1998,  pg.  589). 

Takeoff  Angle  (Degrees): 

Number  of  Bubble  Pulses:  The  number  of  bubble  pulses,  which  can  continue  until  the 
bubble  rises  and  breaks  the  surface. 

Bubble  Pulse  Weighting:  Exponential  time  constants  for  each  bubble  pulse  (?).  Five 
edit  boxes  are  provided,  one  for  each  of  up  to  five  potential  bubble  pulses. 

Save  Parameters:  If  selected  then  the  inputs  in  the  Theoretical  Cepstrum  Parameters 
GUI  are  saved  to  use  in  the  inversion.  Saving  these  parameters  is  optional.  If  not 
selected,  then  the  default  parameters  will  be  used.  Once  the  parameters  are  saved,  they 
will  become  the  new  defaults,  being  saved  as  theo.mat  in  the 
/matlab/toolbox/Cepstrum/MATFULES  directory. 

Close:  If  selected,  closes  Theoretical  Cepstrum  Parameters  window. 
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Figure  A-8.  Theoretical  Cepstrum  Parameters  Window 


A2  EXAMPLE 

The  following  is  a  complete  example  that  includes  processing  stacked  cepstra  from 
selected  data.  The  files  saved  from  processing  the  data  cepstrum  are  then  used  for  the 
inversion.  The  stacked  data  cepstrum  is  compared  to  model  cepstra  to  find  the  best 
match.  This  best  match  provides  the  user  with  the  depth  and  yield  of  the  presumed 
underwater  explosion.  Please  note  that  the  terms  in  bold  are  described  in  Section  A1 . 

A2.1  Location  of  Files. 

The  files  used  for  this  example  are  located  in  the  locations  detailed  below.  To  compute 
the  data  cepstrum,  the  following  files  are  needed: 

*.w  Waveform  File 
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*. affiliation  Network  station  affiliations 

*. arrival  Summary  information  on  a  seismic  arrival 

*. assoc  Data  associating  arrivals  with  origins 

*. origin  Data  on  event  location  and  confidence  bounds 

*.wfdisc  Waveform  file  header  and  descriptive  information 

*.wftag  Waveform  mapping  file 

Note:  In  this  example,  the  .w  file,  and  also  the  additional  files  needed  to  process  the  data 
cepstrum  have  the  prefix  20610545.  The  .w  file  and  the  additional  files  do  not  need  to 
have  the  same  prefix,  although  the  additional  files  all  need  to  be  named  with  a  common 
prefix. 

The  *.w  file  used,  as  well  as  the  additional  files  needed  to  compute  the  cepstrum  are 
located  in  (this  folder  is  arbitrary): 

C  :/CTBT/Data/Dead_Sea/991 1 1  IDeadSea/ 

This  C:/CTBT/Data/Dead_Sea/991 11  IDeadSea/  folder  contains  a  subfolder  called 
/cepstra_files/  It  is  imperative  that  the  folder  which  contains  your  data  contains  a 
subfolder  named  /cepstraJTdes/  in  order  to  save  computed  cepstra. 

C:/CTBT/Data/Dead_Sea/991111DeadSea/cepstra_files/ 

This  folder  is  where  the  calculated  cepstra  files  will  be  saved. 

The  location  of  the  Cepstral  Modeling  and  Inversion  program  is  in  the  matlab  directory: 

C:/matlabR12/toolbox/Cepstrum/ 

This  /Cepstrum  folder  also  contains  a  /MATFILES  folder,  where  the  default  inputs  to  the 
GUI  may  be  saved.  The  /Cepstrum  folder  has  also  been  added  to  the  matlab  path. 
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C:/matlabRl  2/tooIbox/Cepstru  m/M  ATFILES/ 


A2.2  DATA  CEPSTRUM. 

The  data  cepstrum  is  computed  from  inputs  from  the  Stations  and  Channels  GUI  win¬ 
dow,  and  also  the  Process  Cepstrum  GUI  window 

A2.2.1  Stations  and  Channels  —  Selecting  Waveform,  Stations  and  Channels  for 
Computing  Data  Cepstrum 

On  the  matlab  command  line,  the  command  cep_ana  is  typed.  This  opens  the  Modeling 
and  Inversion  GUI  shown  in  Figure  .A-l.  Please  note  that  the  folder  containing  the  cep- 
stral  inversion  algorithm  has  been  added  to  the  matlab  path.  Therefore,  it  can  be  run  from 
any  matlab  directory.  First,  it  is  desired  to  compute  a  cepstrum  from  data.  Therefore,  the 
Parameters  pull  down  menu  is  chosen,  and  the  Stations  and  Channels  option  is  se¬ 
lected,  as  shown  in  Figure  A-9. 

The  Stations  and  Channels  GUT  window  opens,  with  blank  inputs,  as  shown  in  Figure 
A-10.  Click  on  Waveform  Files...  pushbutton.  This  opens  the  Browse  Window, 

Choose  .w  File  to  Analyze  (Figure  A-l  1).  The  Look  In  directory  is  changed  to 
C:/CTBT/Data/Dead  _Sea/991 1 1  IDeadSea/.  The  browse  window  filters  on  .w  files.  The 
desired  2061 0545. w  file  is  selected  with  the  mouse.  Open  may  now  be  pressed  to  open 
the  file,  or  the  file  may  be  double  clicked  to  open. 

The  Stations  and  Channels  GUI  window  gets  filled  in  with  the  waveform  file  name. 

Also  the  Prefix  for  Files  edit  box  is  filled  in  with  the  prefix  of  the  selected  .w  file.  This 
is  the  prefix  for  the  additional  files,  besides  the  .w  file  needed  to  process  the  data.  If  the 
prefix  for  these  additional  files  differs  from  the  prefix  for  the  .w  file,  this  Prefix  for  Files 
edit  box  must  be  changed  to  reflect  the  prefix  for  the  additional  files.  The  Select  All  but¬ 
ton  is  selected.  Therefore  all  of  stations  corresponding  to  the  input  waveform  are  also 
selected.  Clear  All  clears  the  selection  of  the  stations. 

The  program  displays  the  stations  and  channels  associated  with  the  chosen  ,w  file.  Vari- 
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Various  combinations  of  stations  and  channels  may  be  selected  for  processing  the  data 
cepstrum.  However,  at  least  one  station  and  a  corresponding  channel  must  be  selected. 
Once  the  desired  inputs  are  chosen,  push  Save  Selection.  This  saves  the  inputs  for  use  in 
processing  the  data  cepstrum,  and  also  saves  the  input  .w  file  and  location,  to  be  used  as  a 
subsequent  default  as  the  file  file_loc.mat 
inC:/matlabR12/toolbox/Cepstrum/MATFILES. 

For  this  example,  the  bz  channels  of  stations  EEL  and  MRNI  were  selected  for  processing, 
as  shown  Figure  A- 12. 


Figure  A-9.  Modeling  and  Inversion 
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Figure  A- 10.  Initial  Stations  and  Channels  Display 


Figure  A-ll.  Waveform  Files . Browse  Window 
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Figure  A-12.  Selected  Stations  and  Channels 


A2.3  PROCESS  CEPSTRUM  -  VARIABLE  INPUTS  FOR  COMPUTING  DATA 
CEPSTRUM. 

In  the  Modeling  and  Inversion  GUI  window  (Figure  A-l),  from  the  Select  menu, 
choose  Data  Cepstrum,  as  shown  in  Figure  A- 13. 

The  Process  Cepstrum  window  opens,  with  the  *.w  file  chosen  in  the  Stations  and 
Channels  window  displayed  in  the  File  Processing...  display.  The  Inprefix  for  Files 
from  the  Stations  and  Channels  window  is  displayed  in  the  Inprefix  box.  Note  that 
although  browsing  for  different  *.w  files  is  possible  by  clicking  on  the  File  Processing... 
pushbutton,  changes  of  the  file  used  for  processing  must  be  made  in  the  Stations  and 
Channels  window. 
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Defaults  are  used  to  fill  in  the  variable  values  for  the  remaining  inputs.  The  Save  Files 
edit  box  is  filled  in  with  the  path  of  the  *.w  file  from  the  Stations  and  Channels 
window,  with  the  /cepstra_files/  folder  appended.  The  location  of  where  to  save  the 
cepstra  files  may  be  changed.  However,  it  should  be  saved  in  a  folder  by  the  name 
/cepstra_files/.  Figure  A-14  shows  the  Process  Cepstrum  window  with  default  values. 

To  change  values  of  the  inputs  for  processing  the  cepstrum,  simply  change  the  values  of 
the  inputs.  The  values  displayed  on  the  Process  Cepstrum  window  will  be  used  to 
process  the  data  cepstrum.  If  the  altered  inputs  are  to  be  saved  for  use  as  the  future 
defaults,  select  the  pull  down  menu  Utilities,  and  choose  Save  Defaults.  This  will  save 
the  values  to  be  used  as  the  future  defaults  in  the  file  default_data.mat  in  location 
C:/matlabR  1 2/toolbox/ Cepstrum/MATFILES/. 

Figure  A-14  shows  the  Process  Cepstrum  window  with  some  of  the  input  variables 
edited.  To  process  the  cepstrum,  the  pushbutton  Process  Cepstrum  is  pressed. 


Figure  A-13.  Modeling  and  Inversion  --  Data  Cepstrum 


A-21 


J=jn]-Xj 


Utilities 


File  Processing ... 


In  prefix: 


Window  | Hanning  _^_|Length  (s):|  20~ 

Freq  From:j  Ti  To:  |  W~ 


Polynomial 
Fit  Degree 


F  Subtract  Noise 
F  Remove  Instr.  Resp. 

F  Plot  Raw  Spectra 
F  plot  Spectraj 

Cepstrum  Cepstrum  Len.(s) 

C  Stacked  Spectrum 
(*  Stacked  Cepstrum 


1.5 


Save  Files:  j  C:\CTBT\Data\Dead_Sea\ 


Process  Cepstrum 


Close 


Figure  A-14.  Process  Cepstrum  Inputs 


The  inputs  chosen  in  the  Process  Cepstrum  window  produce  Figure  A- 15,  Figure  A-16, 
and  Figure  A-17.  Since  the  Plot  Raw  Spectra  and  Plot  Spectra  options  are  both  chosen, 
the  plots  in  Figure  A- 15  and  Figure  A-16  are  generated.  The  figures  generated  by  the 
cepstral  inversion  program  all  contain  a  heading  consisting  of  the  orid,  date,  time,  stations 
used  to  create  plot  (if  applicable),  and  a  brief  title  describing  the  plot. 

The  individual  spectra,  either  raw  or  processed,  are  plotted  in  a  separate  figure  for  each 
station.  The  spectrum  for  each  individual  phase  for  each  channel  is  plotted  on  the  same 


A-22 


figure  for  each  station.  Since  two  channels  (bn  and  bz)  were  selected  for  both  the  EIL 
and  MRNI  stations  in  the  Process  Cepstrum  window,  two  traces  for  each  phase  are  seen 
plotted  in  the  plots  in  Figure  A-15  and  Figure  A-16.  The  spectra  are  windowed  by  a 
Hanning  window,  and  are  a  length  of  20  seconds,  set  in  Window  and  Length  (s)  in  the 
Process  Cepstrum  window.  The  processed  cepstrum  reflects  the  Freq  From:  and  To: 
values  of  1 .5  to  18  (Hz),  and  also  the  polynomial  degree  fit  subtraction  of  2  from  the 
Polynomial  Fit  Degree  input.  Since  Subtract  Noise  and  Remove  Instr.  Resp.  are  both 
checked,  the  noise  was  subtracted  from  each  of  the  spectrum  phases,  and  the  instrument 
response  was  removed.  The  time  length  of  the  cepstrum  in  Figure  A- 17  reflects  the  user 
input  Cepstrum  Len.(s)  of  1.5  (s).  The  cepstrum  was  computed  by  stacking  all  of  the 
phase  cepstra. 

The  spectra  plots  of  Figure  A-15  and  Figure  A-16  display  the  spectra  of  all  phases 
available  for  the  station.  However,  for  computing  the  stacked  cepstrum  in  Figure  A-17, 
only  the  phases  that  are  available  across  all  of  the  chosen  stations  are  used.  Since  only 
the  Pn  phase  is  available  for  stations  EEL  and  MRNI,  only  this  station  is  used  to  compute 
the  stacked  cepstrum. 

The  cepstra  files  are  saved  in  the 

C:/CTBT/Data/Dead_Sea/991 1 1 1  DeadSea/cepstra_files/  folder  specified  in  the  Save 
Files  edit  box  in  the  Process  Cepstrum  window.  The  saved  files  are: 

2061 0545. mat  Contains  the  variables  used  for  processing  that  are  also  needed  for 
inversion 

20610545. ST. cp  Stacked  Cepstrum 
20610545. Ns. cp  Stacked  Noise  Cepstrum 
20610545. Lg.cp  Stacked  Phase  Cepstrum  (Lg  phase) 
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Figure  A-15.  Raw  and  Processed  Spectra  for  Each  Phase  for  EIL  Station 
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Figure  A- 16.  Raw  and  Processed  Spectra  for  Each  Phase  for  MRNI  Station 
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Figure  A-17.  Cepstra  Computed  with  User  Inputs 


A2.4  CALCULATING  DEPTH  AND  YIELD  BY  INVERSION. 

The  Inversion  of  the  data  cepstrum  to  the  theoretical  cepstrum  is  computed  using  the 
selected  saved  cepstra  files,  and  the  inputs  to  the  Inversion  and  Theoretical  Cepstrum 
Parameters  window.  Inversion  of  the  data  cepstrum  to  theoretical  cepstra  determines 
the  best  fit  of  depth  and  yield  of  the  explosion 

A2.4.1  Inversion 

In  the  Modeling  and  Inversion  window  (Figure  A-l),  from  the  Select  menu,  choose  In¬ 
version,  as  shown  in  Figure  A- 18. 
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The  Inversion  window  opens,  with  blank  inputs  for  the  Saved  Files...  and  the  Inprefix:. 
The  inputs  for  the  variables  filling  in  the  rest  of  the  window  are  defaults  (Figure  A-19). 
Click  on  the  Saved  Files...  pushbutton.  This  opens  the  Browse  Window,  Choose  ,cp  or 
.sp  File  to  Analyze  (Figure  A-20).  The  Look  In  directory  is  changed  to 
C:/CTBT/Data/Dead_Sea/991 1 1  lDeadSea/cepstra_files/.  The  browse  window  filters  on 
.ST.cp  files  (Figure  A-20).  These  files  are  the  saved  stack  cepstra  files  saved  from  the 
Process  Cepstrum  window.  The  desired  206 10545. ST.cp  file  is  selected  with  the 
mouse.  Open  may  now  be  pressed  to  open  the  file,  or  the  file  may  be  double  clicked  to 
open. 

The  Saved  Files...  box  is  filled  in  with  the  .ST.cp  file  chosen.  The  Inprefix  edit  box  is 
filled  in  with  the  prefix  of  the  chosen  cepstra  file.  To  change  values  of  the  parameter 
inputs  for  the  inversion,  simply  change  the  values  of  the  inputs.  The  values  displayed  in 
the  Inversion  window  will  be  used  to  process  the  inversion.  If  the  altered  inputs  are  to 
be  saved  for  use  as  the  future  defaults,  select  the  pull  down  menu  Utilities,  and  choose 
Save  Defaults.  This  will  save  the  values  to  be  used  as  the  future  defaults  in  the  file 
inv.mat  in  location  C:/matlabR12/toolbox/Cepstrum/MATFILES/. 

In  order  to  determine,  and  also  change  the  default  theoretical  parameters,  the  Select 
Theoretical  Parameters  pushbutton  is  pressed  in  the  Inversion  window.  This  opens  the 
Theoretical  Cepstrum  Parameters  window,  with  the  default  parameters  as  inputs 
(Figure  A-22).  Any  of  the  input  parameters  may  be  edited.  To  use  the  changed 
parameters  in  the  current  inversion  calculation,  check  the  Save  Parameters  checkbox  in 
the  Theoretical  Parameters  window.  To  save  the  parameters  for  future  default  settings, 
in  the  Inversion  window,  select  the  pull  down  menu  Utilities,  and  choose  Save  Defaults. 
Note  that  this  will  also  save  the  inputs  in  the  Inversion  window  to  use  as  future  defaults. 
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Figure  A-18.  Modeling  and  Inversion  Program  -  Inversion 
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Figure  A-19.  Initial  Inversion  Window 
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Figure  A-20.  Files  Saved...  Browse  Window 
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Figure  A-21 .  Inversion  Inputs 
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Figure  A-22.  Theoretical  Cepstrum  Parameters  Window  with  Defaults 


A2.4.2  Inversion  using  LI  Norm  Comparison  Method,  and  Range  of  Values 
Optimization  Method 

The  inputs  chosen  in  the  Inversion  window  produce  the  related  figures  (Figure  A-23, 
Figure  A- 24,  Figure  A-25,  Figure  A-26).  First  displayed  is  an  image  of  the  saved  cepstra 
phase  and  stacked  cepstra,  used  in  the  inversion  algorithm  (Figure  A-23).  The  input  of 
the  Saved  Files...  and  Inprefix  boxes  determine  which  saved  cepstra  files  to  use  in  the 
inversion  and  display.  Parameters  input  in  the  Inversion  window  determine  the 
comparison  method  of  the  data  cepstrum  to  the  theoretical  cepstra.  In  this  example,  LI 
Norm  was  chosen  as  the  Comparison  method,  with  a  Range  of  Values  Optimization 
method.  The  Yield(kg):  and  To:  values  of  2000  and  6000,  respectively,  is  the  range  of 
yields  for  the  theoretical  cepstra,  with  a  Step:  of  100. 
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The  data  cepstrum  is  compared  to  theoretical  cepstra,  computed  with  the  inputs  of  the 
Inversion  and  Theoretical  Cepstrum  Parameters  inputs.  Comparison  of  the  data 
cepstrum  with  the  range  of  values  input  into  creating  various  theoretical  cepstra  produces 
corresponding  1  -LI  Norm  values.  The  range  of  values  chosen  in  the  Inversion  window 
is  reflected  in  the  3D  plot  of  1  -LI  Norm  of  the  data  cepstrum  and  theoretical  cepstra  over 
the  input  range  of  values  (Figure  A-24).  A  contour  plot  of  the  range  of  values  1 -LI  Norm 
is  also  created,  shown  in  Figure  A-25.  A  search  is  made  for  the  best  match  of  1-L1  Norm 
between  the  data  cepstrum  and  the  theoretical  cepstra.  The  best  matching  theoretical 
cepstrum  is  plotted,  in  red  dashes,  along  with  the  data  cepstrum  (Figure  A-26).  The 
depth,  yield,  and  match  values  are  listed  in  the  plot  title.  The  second  plot  in  Figure  A-26 
shows  the  plot  of  1 -LI  Norm,  with  depth  constrained  to  the  depth  of  maximizing  the 
comparison.  The  plots  in  Figure  A-26  are  also  shown  for  the  second  highest  match 
between  the  theoretical  and  data  cepstrum. 
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Figure  A-23.  Image  of  Saved  Phase  and  Stack  Cepstra  Used  for  Inversion 
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Figure  A-24.  3D  Range  of  Values 
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Figure  A-25.  Contour  Plot  of  1  -LI  Norm  of  Data  Cepstrum  and  Theoretical  Cepstra 
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Figure  A-26.  Comparison  of  Theoretical  and  Data  Cepstrum,  and  Plot  of  1  -LI  Norm  vs. 
Yield 

A2.4.3  Comparison  Methods 

The  example  in  section  A2.8  used  the  LI  Norm  Comparison  method  for  finding  the 
depth  and  yield  of  the  explosion.  Cross  Corr,  or  L2  Norm  may  also  be  used  as  the 
Comparison  method  between  the  data  cepstrum  and  various  theoretical  cepstra.  The 
same  types  of  figures  as  those  described  above  will  be  created,  and  labeled  with  the 
appropriate  Comparison  choice. 

A2.4.4  Optimization  Methods 

The  example  in  A2.8  used  the  Range  of  Values  Optimization  method.  All 
combinations  between  Comparison  and  Optimization  can  be  used  to  find  the  depth  and 
yield  of  the  explosion. 

Choosing  Downhill  Simplex  as  the  Optimization  method  changes  the  Inversion 
window  to  look  as  that  in  Figure  A- 27.  The  input  required  for  Yield(kg):  is  a  single 
value,  which  is  the  starting  point  of  the  Downhill  Simplex  algorithm.  The  single 
Depth(m):  value  also  is  a  starting  point  for  the  algorithm.  In  this  example,  the  default 
yield  and  depth  were  changed  to  match  the  best  match  yield  and  depth  determined  by  the 
inversion  algorithm  for  the  range  of  values  example  in  A2.4.2. 
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As  for  each  inversion,  the  first  figure  created  displays  the  data  cepstra,  as  in  Figure  A-23. 
Each  value  of  depth,  yield,  and  surface  reflection  coefficient  the  downhill  simplex 
algorithm  chose,  is  marked  with  a  blue  ‘x’  in  Figure  A-28.  Also  shown  is  the  best  match, 
determined  by  the  downhill  simplex  method  for  the  comparison  method  of  choice 
between  the  data  and  theoretical  cepstra. 

Since  Plot  Theoretical  Model  was  checked,  Figure  A-29  was  created.  This  shows  the 
components  of  the  model  theoretical  cepstrum  for  the  depth  and  yield  values  of  the  best 
matching  cepstrum. 

Choosing  the  Simulated  Annealing  Optimization  method  uses  an  algorithm  that 
attempts  to  find  a  global  maximum  match  between  the  data  cepstrum  and  theoretical 
cepstra.  When  this  Optimization  method  is  chosen,  the  Temps:  and  Trials:  edit  boxes 
become  active.  These  may  be  changed  to  reflect  the  desired  number  of  trial  and 
temperatures  for  the  algorithm.  In  this  example  Depth  and  Yield  from  Data  Cepstrum 
was  selected.  This  option  allows  the  user  to  click  on  the  bottom  reflection  and  first 
bubble  pulse  in  the  data  cepstrum  figure  (Figure  A-23),  which  is  plotted  for  each 
inversion  computed.  The  user  places  the  cursor  over  the  desired  location,  a  cross 
designates  the  locadon,  and  the  user  clicks  (Figure  A-30). 

The  results  of  the  simulated  annealing  algorithm  search  are  plotted,  and  the  best  depth 
and  yield  found. 
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Figure  A-27.  Inversion  with  Optimization  Method  Downhill  Simplex 
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Figure  A-28.  Downhill  Simplex  Best  Value  Search,  and  Best  Match  of  Data  and 
Theoretical  Cep  strum 
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Figure  A-29.  Theoretical  Model 
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Figure  A-30.  Determining  Starting  Depth  and  Yield  from  User  Inputs 
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Figure  A-31 .  Simulated  Annealing,  and  Best  Match  of  Data  and  Theoretical  Cepstra 
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1408  NE  45TH  ST  0201 
SEATTLE,  WA  98105 

ATTN:  DR.  T.  AHERN 

MASSACHUSETTS  INSTITUTE  OF  TECHNOLOGY 
EARTH  RESOURCES  LABORATORY 
42  CARLETON  STREET 
CAMBRIDGE,  MA  02142 

ATTN:  DR.  W.  RODI 
ATTN:  PROF.  M.N.  TOKSOZ 

MICHIGAN  STATE  UNIVERSITY  LIBRARY 
450  ADMINISTRATION  BUILDING 
EAST  LANSING,  Ml  48824 
ATTN:  K.  FUJITA 

NEW  MEXICO  STATE  UNIVERSITY 
DEPARTMENT  OF  PHYSICS 
LAS  CRUCES,  NM  88003 

ATTN:  PROF.  J.  Nl 
ATTN:  PROF.  T.  HEARN 

NORTHWESTERN  UNIVERSITY 
DEPARTMENT  OF  GEOLOGICAL  SCIENCES 
1847  SHERIDAN  RD 
EVANSTON,  IL  60208 

ATTN:  PROF.  E.  OKAL 

PENNSYLVANIA  STATE  UNIVERSITY 
GEOSCIENCES  DEPARTMENT 
403  DEIKE  BUILDING 
UNIVERSITY  PARK,  PA  16802 

ATTN:  PROF.  C.  AMMON 
ATTN:  PROF.  S.  ALEXANDER 
ATTN:  DR.  A.  NYBLADE 

SAN  DIEGO  STATE  UNIVERSITY 
DEPARTMENT  OF  GEOLOGICAL  SCIENCES 
SAN  DIEGO,  CA  92182 

ATTN:  PROF.  S.  M.  DAY 

SOUTHERN  METHODIST  UNIVERSITY 
DEPARTMENT  OF  GEOLOGICAL  SCIENCES 
P.O.  BOX  750395 
DALLAS,  TX  75275 

ATTN:  B.  STUMP 
ATTN:  E.  HERRIN 
ATTN:  P.  GOLDEN 

UNIVERSITY  OF  HAWAII-  MANOA 
P.O.  BOX  1599 

KAILUA-KONA,  HI  96745  1599 

ATTN:  DR.  M.  A.  GARCES 

UNIVERSITY  OF  MISSISSIPPI 
1  COLISEUM  DRIVE 
UNIVERSITY,  MS  38677 

ATTN:  PROF.  H.  BASS 
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UNIVERSITY  OF  SOUTHERN  CALIFORNIA 
520  SEAVER  SCIENCE  CENTER 
UNIVERSITY  PARK 
LOS  ANGELES,  CA  90089  0483 

ATTN:  PROF.  C.  G.  SAMMIS 
ATTN:  PROF.  T.  JORDAN 

UNIVERSITY  OF  WISCONSIN  MADISON 
1215  W  DAYTON  ST 
MADSION,  Wl  53706  1600 

ATTN:  DR.  C.  THURBER 

ST  LOUIS  UNIVERSITY 
EARTH  &  ATMOSPHERIC  SCIENCES 
STATION  3507  LACLEDE  AVE 
ST  LOUIS,  MO  63103 

ATTN:  PROF.  B.  J.  MITCHELL 
ATTN:  PROF.  R.  HERRMAN 

UNIVERSITY  OF  MEMPHIS 
3890  CENTRAL  AVE 
MEMPHIS,  TN  38152 

ATTN:  DR.  J.  PUJOL 

UNIVERSITY  OF  MEMPHIS 
3876  CENTRAL  AVE 
MEMPHIS,  TN  38152 

ATTN:  DR.  C.  LANGSTON 

UNIVERSITY  OF  TEXAS  AUSTIN 
IGS  130 

AUSTIN,  TX  78712 

ATTN:  DR.  J.  PULLIAM 

UNIVERSITY  OF  TEXAS  AUSTIN 
IGS  131 

AUSTIN,  TX  78712 

ATTN:  DR.  M.  SEN 

UNIVERSITY  OF  TEXAS  EL  PASO 
DEPT  OF  GEOLOGICAL  SCIENCES 
EL  PASO,  TX  79901 

ATTN:  PROF.  G.  KELLER 
ATTN:  DR.  D.  DOSER 
ATTN:  DR.  A.  VELASCO 
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AUSTRALIAN  GEOLOGICAL  SURVEY 
ORGANIZATION 

CORNER  OF  JERRAGOMRRRA  & 
N1NDMARSH  DRIVE 
CANBERRA,  ACT  2609 
AUSTRALIA 

ATTN:  D.  JEPSON 

GEOPHYSICAL  INSTITUTE  OF  ISRAEL 
POB  182 

LOD,  7100  ISRAEL 

ATTN:  DR.  Y.  GITTERMAN 
ATTN:  DR.  A.  SHAPIRA 

GEOLOGICAL  SURVEY  OF  CANADA 
7  OBSERVATORY  CRESCENT 
OTTAWA  K1A  0Y3  ONT 
CANADA 

ATTN:  C.  WOODGOLD 

I.R.I.G.M.  B.P.  68 

38402  ST.  MARTIN  D'HERES 

CEDEX,  FRANCE 

ATTN:  DR.  M.  BOUCHON 

MINISTRY  OF  DEFENSE 
PROCUREMENT  EXEUTIVE 
BLACKNESS,  BRIMPTON 
READING  FG7-4RS  ENGLAND 

ATTN:  DR.  D.  BOWERS 

NTNF/NORSAR 
P.O.  BOX  51 

N-2007  KJELLER,  NORWAY 

ATTN:  DR.  F.  RINGDAL 
ATTN:  T.  KVAERNA 
ATTN:  S.  MYKKELTVEIT 
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